!/===========================================================================/
! Copyright (c) 2007, The University of Massachusetts Dartmouth 
! Produced at the School of Marine Science & Technology 
! Marine Ecosystem Dynamics Modeling group
! All rights reserved.
!
! FVCOM has been developed by the joint UMASSD-WHOI research team. For 
! details of authorship and attribution of credit please see the FVCOM
! technical manual or contact the MEDM group.
!
! 
! This file is part of FVCOM. For details, see http://fvcom.smast.umassd.edu 
! The full copyright notice is contained in the file COPYRIGHT located in the 
! root directory of the FVCOM code. This original header must be maintained
! in all distributed versions.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
! AND ANY EXPRESS OR  IMPLIED WARRANTIES, INCLUDING,  BUT NOT  LIMITED TO,
! THE IMPLIED WARRANTIES OF MERCHANTABILITY AND  FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED.  
!
!/---------------------------------------------------------------------------/
! CVS VERSION INFORMATION
! $Id$
! $Name$
! $Revision$
!/===========================================================================/

MODULE RRKVAL
# if defined (RRKF)
   
   USE CONTROL
   IMPLICIT NONE
  
   REAL(DP),ALLOCATABLE :: STTEMP0(:)
   REAL(DP),ALLOCATABLE :: STTEMP1(:)
   REAL(DP),ALLOCATABLE :: STEOF(:)
   REAL(DP),ALLOCATABLE :: SDEOF(:)
   REAL(DP),ALLOCATABLE :: TRANS(:,:)
   REAL(SP) ::  OBSERR_EL         !!EL OBSERVATION ERROR SPECIFIED
   REAL(SP) ::  OBSERR_UV         !!UV OBSERVATION ERROR SPECIFIED
   REAL(SP) ::  OBSERR_T          !!TEMPERATURE OBSERVATION ERROR SPECIFIED
   REAL(SP) ::  OBSERR_S          !!SALINITY OBSERVATION ERROR SPECIFIED   
   REAL(SP),ALLOCATABLE :: RRKEL(:)
   REAL(SP),ALLOCATABLE :: RRKU(:,:)
   REAL(SP),ALLOCATABLE :: RRKV(:,:)
   REAL(SP),ALLOCATABLE :: RRKT(:,:)
   REAL(SP),ALLOCATABLE :: RRKS(:,:)

   REAL(SP),ALLOCATABLE :: RRKEL3(:)
   REAL(SP),ALLOCATABLE :: RRKU3(:,:)
   REAL(SP),ALLOCATABLE :: RRKV3(:,:)
   REAL(SP),ALLOCATABLE :: RRKT3(:,:)
   REAL(SP),ALLOCATABLE :: RRKS3(:,:)

   REAL(SP), ALLOCATABLE :: RRKU_REF(:,:,:)
   REAL(SP), ALLOCATABLE :: RRKV_REF(:,:,:)
   REAL(SP), ALLOCATABLE :: RRKEL_REF(:,:)
   REAL(SP), ALLOCATABLE :: RRKT_REF(:,:,:)
   REAL(SP), ALLOCATABLE :: RRKS_REF(:,:,:)
   
   REAL(DP) ::  ELSD, USD, TSD, SSD
   REAL(DP),ALLOCATABLE :: STSD(:)
   REAL(DP),ALLOCATABLE :: SEOF(:,:)
   REAL(DP),ALLOCATABLE :: STMEAN(:)
   INTEGER  nrange
   INTEGER  ::  NLOC=0              !!Number of observation locations
   LOGICAL  ::  EL_ASSIM          !!OPTION FOR CHOSING ELEVATION AS ASSIMILATION VARIABLES
   LOGICAL  ::  UV_ASSIM          !!OPTION FOR CHOSING CURRENT AS ASSIMILATION VARIABLES
   LOGICAL  ::  T_ASSIM           !!OPTION FOR CHOSING TEMPERATURE AS ASSIMILATION VARIABLES
   LOGICAL  ::  S_ASSIM           !!OPTION FOR CHOSING SALINITY AS ASSIMILATION VARIABLES

   LOGICAL  ::  EL_OBS            !!OPTION FOR ELEVATION OBSERVATION DATA
   LOGICAL  ::  UV_OBS            !!OPTION FOR CURRENT OBSERVATION DATA
   LOGICAL  ::  T_OBS             !!OPTION FOR TEMPERATURE OBSERVATION DATA
   LOGICAL  ::  S_OBS             !!OPTION FOR SALINITY OBSERVATION DATA
# endif
END MODULE RRKVAL

MODULE MOD_RRK
   USE MOD_INPUT, only : NC_START
   USE RRKVAL
   USE CONTROL
   USE MOD_UTILS
   USE MOD_NCTOOLS
   IMPLICIT NONE
   SAVE
   
   INTEGER  ::  RRK_EOFCONTR
   CHARACTER(LEN=80) :: REF_START_DATE
   CHARACTER(LEN=80) :: REF_END_DATE
   CHARACTER(LEN=80) :: RRK_START_DATE
   CHARACTER(LEN=80) :: RRK_END_DATE

   TYPE(TIME) :: REF_START_TIME
   TYPE(TIME) :: REF_END_TIME
   TYPE(TIME) :: RRK_START_TIME
   TYPE(TIME) :: RRK_END_TIME

   TYPE(TIME) :: RRK_CYC

   CHARACTER(LEN=80) :: RRK_ASSIM_INTERVAL
   TYPE(TIME) :: RRK_INTERVAL


   INTEGER      REF_INT           !!GLOBAL NUMBER OF THE READING FILE INTERVALS 
   INTEGER      RRK_NOBSMAX
   INTEGER      RRK_OPTION        !!OPTION 1 FOR BAROTROPIC CASE; OPTION 2 FOR BAROCLINIC CASE
   INTEGER      RRK_NEOF          !!NUMBER OF THE EOF  
   REAL(SP) ::  RRK_PSIZE         !!PERTURBATION SIZE  
   REAL(SP) ::  RRK_PSCALE        !!PSEUDO MODEL ERROR    
   REAL(SP) ::  RRK_RSCALE        !!SCALE FACTOR APPLIED TO ONE STANDARD DEVIATION FOR R

   LOGICAL  :: RRK_ON

# if defined (RRKF)
   INTEGER  ::  ICYC
!   LOGICAL  ::  EL_ASSIM          !!OPTION FOR CHOSING ELEVATION AS ASSIMILATION VARIABLES
!   LOGICAL  ::  UV_ASSIM          !!OPTION FOR CHOSING CURRENT AS ASSIMILATION VARIABLES  
!   LOGICAL  ::  T_ASSIM           !!OPTION FOR CHOSING TEMPERATURE AS ASSIMILATION VARIABLES
!   LOGICAL  ::  S_ASSIM           !!OPTION FOR CHOSING SALINITY AS ASSIMILATION VARIABLES
!
!   LOGICAL  ::  EL_OBS            !!OPTION FOR ELEVATION OBSERVATION DATA
!   LOGICAL  ::  UV_OBS            !!OPTION FOR CURRENT OBSERVATION DATA  
!   LOGICAL  ::  T_OBS             !!OPTION FOR TEMPERATURE OBSERVATION DATA
!   LOGICAL  ::  S_OBS             !!OPTION FOR SALINITY OBSERVATION DATA

   INTEGER      INORRK            !!FILE I/O PIPE NUMBER 


   TYPE(NCFILE), POINTER :: NCF_RRKFDATA

   TYPE(NCFILE), POINTER :: NCF_RRKRE
   TYPE(NCFILE), POINTER :: NCF_RRKFCT
   LOGICAL :: LOCAL_DISK
   CHARACTER(LEN=80), parameter :: Local_file ="/scratch/rrk_restart.nc"
   CHARACTER(LEN=80), parameter :: group_file ="rrktemp/rrk_restart.nc"

   CHARACTER(LEN=80), parameter :: Local_file1 ="/scratch/rrk_forecast.nc"
   CHARACTER(LEN=80), parameter :: group_file1 ="rrktemp/rrk_forecast.nc"

   INTEGER(ITIME) :: RKINT_COUNT, RKINT_START, RKINT_END


   INTERFACE READRESTART
      MODULE PROCEDURE READRESTART_TIME
      MODULE PROCEDURE READRESTART_STATE
   END INTERFACE

   PRIVATE READRESTART_TIME
   PRIVATE READRESTART_STATE


   NAMELIST /NML_RRKF/        &
         & RRK_ON,            &
         & REF_START_DATE,    &
         & REF_END_DATE,      &
         & RRK_START_DATE,    &
         & RRK_END_DATE,      &
         & RRK_ASSIM_INTERVAL,&
         & RRK_NOBSMAX,       &
         & RRK_NEOF,          &
         & RRK_PSIZE,         &
         & RRK_PSCALE,        &
         & RRK_RSCALE,        &
         & EL_ASSIM,          &
         & EL_OBS,            &
         & UV_ASSIM,          &
         & UV_OBS,            &
         & T_ASSIM,           &
         & T_OBS,             &
         & S_ASSIM,           &
         & S_OBS,             &
         & LOCAL_DISK

   CONTAINS
   
     SUBROUTINE NAME_LIST_INITIALIZE_RRKF
       IMPLICIT NONE
       REF_START_DATE = "RRK REFERENCE START AND END TIME FOR EOF CALCULATION"
       REF_END_DATE   = "Date and Time are specified as a string (example '2007-11-05 00:00:00')"  
       RRK_START_DATE = "RRK ASSIMILATION START AND END TIME"
       RRK_END_DATE   = "For an idealized case specify 'seconds=(flt)','days=(flt)', or 'cycles=(int)'" 
       RRK_ASSIM_INTERVAL= "A length of time: 'seconds= ','days= ', or 'cycles= '"   
       RRK_NOBSMAX    = -1  
       RRK_NEOF       = -1  
       RRK_PSIZE      = -1  
       RRK_PSCALE     = -1.0_SP  
       RRK_RSCALE     = -1.0_SP  
       RRK_ON         = .FALSE.
       EL_ASSIM       = .FALSE.  
       EL_OBS         = .FALSE.  
       UV_ASSIM       = .FALSE.  
       UV_OBS         = .FALSE. 
       T_ASSIM        = .FALSE.   
       T_OBS          = .FALSE.  
       S_ASSIM        = .FALSE.  
       S_OBS          = .FALSE.
       LOCAL_DISK     = .FALSE.

     END SUBROUTINE NAME_LIST_INITIALIZE_RRKF
     
     SUBROUTINE NAME_LIST_PRINT_RRKF
       IMPLICIT NONE

       write(UNIT=IPT,NML=NML_RRKF)

     RETURN
     END SUBROUTINE NAME_LIST_PRINT_RRKF
     
     SUBROUTINE NAME_LIST_READ_RRKF
       IMPLICIT NONE
       CHARACTER(LEN=120) :: FNAME
       INTEGER :: IOS

       FNAME = TRIM(INPUT_DIR)//"/"//trim(casename)//"_assim_rrkf.nml"
       
       CALL FOPEN(KFUNIT,trim(FNAME),'cfr')

       ! Read RRKF parameters
       READ(UNIT=KFUNIT, NML=NML_RRKF,IOSTAT=ios)
       if(ios .NE. 0 ) then
          if(DBG_SET(dbg_log)) write(UNIT=IPT,NML=NML_RRKF)
          Call Fatal_Error("Can Not Read NameList NML_RRKF from file: "//trim(FNAME))
       end if
       
       REWIND(KFUNIT)
       
       if(DBG_SET(dbg_scl)) &
            & write(IPT,*) "Read_Name_List: NML_RRKF"
       
       if(DBG_SET(dbg_scl)) &
            & write(UNIT=IPT,NML=NML_RRKF)
       
       CLOSE(NMLUNIT)
       
       RETURN
     END SUBROUTINE NAME_LIST_READ_RRKF
     
     SUBROUTINE RRK_SET_TIME
       USE MOD_SET_TIME
       IMPLICIT NONE
       CHARACTER(LEN=4) :: FLAG,BFLAG,EFLAG
       INTEGER(ITIME):: begins, ends,steps
       integer :: status
# if defined (DATA_ASSIM)
   FVCOM_RUN_MODE = FVCOM_RRKF_WITH_SSA
# else
   FVCOM_RUN_MODE = FVCOM_RRKF_WITHOUT_SSA
# endif

       
       ! GET THE RRK ASSIMILATION INTERVAL PERIOD
       CALL IDEAL_TIME_STRING2TIME(RRK_ASSIM_INTERVAL,FLAG,RRK_INTERVAL,steps)
       IF (FLAG == 'time') THEN ! IF START AND END TIME WERE SPECIFIED

          IF(MOD(RRK_INTERVAL,IMDTI) /= ZeroTime) THEN
             CALL FATAL_ERROR("RRK_ASSIM_INTERVAL must be an integer number of internal time steps!")
          END IF
          
          steps = seconds(RRK_INTERVAL)/seconds(IMDTI)


       ELSEIF(FLAG == 'step') THEN ! IF SPECIFIED AS A NUMBER OF STEPS
          
          RRK_INTERVAL = steps * IMDTI
          
       END IF

       RKInt_start=0
       RKint_end = steps


       SELECT CASE(USE_REAL_WORLD_TIME)
       CASE(.TRUE.)

          REF_START_TIME = READ_DATETIME(REF_START_DATE,DATE_FORMAT,TIMEZONE,status)
!!$          if (.not. status) &
          if (status == 0) &
               & Call Fatal_Error("Could not read the date string REF_START_DATE: "//trim(REF_START_DATE))


          ! GET THE END TIME
          REF_END_Time = READ_DATETIME(REF_END_DATE,DATE_FORMAT,TIMEZONE,status)
!!$          if (.not. status) &
          if (status == 0) &
               & Call Fatal_Error("Could not read the date string REF_END_DATE: "&
               &//trim(REF_END_DATE))

          ! SANITY CHECK
          if(REF_Start_Time .GE. REF_End_Time) &
               & Call Fatal_Error("RRKF REF_Start_Date exceeds or equal to REF_End_Date")
          
          RRK_START_TIME = READ_DATETIME(RRK_START_DATE,DATE_FORMAT,TIMEZONE,status)
!!$          if (.not. status) &
          if (status == 0) &
               & Call Fatal_Error("Could not read the date string RRK_START_DATE: "//trim(RRK_START_DATE))
          
          
          ! GET THE END TIME
          RRK_END_Time = READ_DATETIME(RRK_END_DATE,DATE_FORMAT,TIMEZONE,status)
!!$          if (.not. status) &
          if (status == 0) &
               & Call Fatal_Error("Could not read the date string RRK_END_DATE: "&
               &//trim(RRK_END_DATE))

          ! SANITY CHECK
          if(RRK_Start_Time .GE. RRK_End_Time) &
               & Call Fatal_Error("RRKF RRK_Start_Date exceeds or equal to RRK_End_Date")
          


       CASE (.FALSE.)

          ! GET THE REFERENCE START AND END INFORMATION
          CALL IDEAL_TIME_STRING2TIME(REF_START_DATE,BFLAG,REF_START_TIME,begins)

          CALL IDEAL_TIME_STRING2TIME(REF_END_DATE,EFLAG,REF_End_TIME,ends)

          ! SANITY CHECK
          IF (BFLAG /= EFLAG) CALL FATAL_ERROR&
               ('IDEALIZED MODEL TIME SPECIFICATION IS INCORRENT',&
               &'RRK REF BEGIN AND END CAN BE IN EITHER CYCLES OR TIME BUT NOT MIXED',&
               & trim(REF_start_date),trim(ref_end_date) )
          
          IF (EFLAG == 'time') THEN ! IF START AND END TIME WERE SPECIFIED

             
             ! DO NOTHING

          ELSE IF(EFLAG == 'step') THEN ! IF START AND END IINT WERE SPECIFIED


             REF_START_TIME = StartTime + IMDTI*(begins - ISTART)

             REF_END_TIME = StartTime + IMDTI*(ends - ISTART)

             
          ELSE
             CALL FATAL_ERROR&
                  &('IDEAL_TIME_STRING2TIME returned invalid flag for rrk ref date and time')

          END IF

          ! GET THE RRK START AND END INFORMATION
          CALL IDEAL_TIME_STRING2TIME(RRK_START_DATE,BFLAG,RRK_START_TIME,begins)

          CALL IDEAL_TIME_STRING2TIME(RRK_END_DATE,EFLAG,RRK_End_TIME,ends)

          ! SANITY CHECK
          IF (BFLAG /= EFLAG) CALL FATAL_ERROR&
               ('IDEALIZED MODEL TIME SPECIFICATION IS INCORRENT',&
               &'RRK BEGIN AND END CAN BE IN EITHER CYCLES OR TIME BUT NOT MIXED',&
               & trim(ref_start_date),trim(ref_end_date) )
          
          IF (EFLAG == 'time') THEN ! IF START AND END TIME WERE SPECIFIED

             
             ! DO NOTHING

          ELSE IF(EFLAG == 'step') THEN ! IF START AND END IINT WERE SPECIFIED


             RRK_START_TIME = StartTime + IMDTI*(begins - ISTART)

             RRK_END_TIME = StartTime + IMDTI*(ends - ISTART)

             
          ELSE
             CALL FATAL_ERROR&
                  &('IDEAL_TIME_STRING2TIME returned invalid flag for rrk ref date and time')

          END IF



       END SELECT


       IF(MOD(REF_START_TIME-StartTime,IMDTI) /= ZeroTime) THEN
          CALL FATAL_ERROR("REF_START_TIME must be an integer number of internal time steps from the start time!")
       END IF
       
       IF(MOD(REF_END_TIME-StartTime,IMDTI) /= ZeroTime) THEN
          CALL FATAL_ERROR("REF_END_TIME must be an integer number of internal time steps from the start time!")
       END IF
       
       
       IF(MOD(RRK_START_TIME-StartTime,IMDTI) /= ZeroTime) THEN
          CALL FATAL_ERROR("RRK_START_TIME must be an integer number of internal time steps from the start time!")
       END IF
       
       IF(MOD(RRK_END_TIME-StartTime,IMDTI) /= ZeroTime) THEN
          CALL FATAL_ERROR("RRK_END_TIME must be an integer number of internal time steps from the start time!")
       END IF
       
       
     END SUBROUTINE RRK_SET_TIME


     SUBROUTINE RRK_SET_STARTUP
       USE MOD_SET_TIME
       USE MOD_INPUT
       IMPLICIT NONE
       character(len=160) :: pathnfile
       TYPE(NCFILE), POINTER :: NCF

       OPEN(75,FILE='./err.dat')

       !RESTART FILE
       RESTART_FILE_NAME = trim(casename)//"_restart_0001.nc"
       pathnfile= trim(OUTPUT_DIR)//TRIM(RESTART_FILE_NAME)
       
       CALL SEARCH_FOR_LAST_MATCHING_NAME(PATHNFILE)
       
       ! OPEN THE FILE AND LOAD FOR STARTUP
       NCF => NEW_FILE()
       NCF%FNAME=trim(pathnfile)
       Call NC_OPEN(NCF)
       CALL NC_LOAD(NCF)
       NC_START => NCF
       
       Nullify(NCF)
       CALL SET_STARTUP_FILE_STACK(REF_START_TIME)

!       STARTUP_TYPE = STARTUP_TYPE_CRASHRESTART

!       STARTUP_UV_TYPE = STARTUP_TYPE_SETVALUES

!       STARTUP_TURB_TYPE = STARTUP_TYPE_SETVALUES

!       STARTUP_TS_TYPE = STARTUP_TYPE_SETVALUES

     END SUBROUTINE RRK_SET_STARTUP


!!=====================================================================================
!!
!! READ IN THE SIMULATION DATA AND CALCULATE THE EOF
!!
!!=====================================================================================

   
   SUBROUTINE RRK_REF 
   USE CONTROL
   USE RRKVAL
   USE MOD_NCTOOLS
   IMPLICIT NONE

   TYPE(NCFILE), POINTER :: NCF
   TYPE(TIME) :: tTEST
   LOGICAL :: FOUND
   INTEGER :: NTIMES
   
   INTEGER IDUMMY, I, K, J
   INTEGER SS_DIM
   INTEGER STDIM
!   REAL(DP) ::  ELSD, USD, TSD, SSD
   INTEGER LWORK4, LDVT
   REAL(DP),ALLOCATABLE :: WORK4(:)
   REAL(DP) :: VT
   REAL(DP),ALLOCATABLE :: RKSF(:,:)
   REAL(DP),ALLOCATABLE :: RKSF1(:,:)
!   REAL(DP),ALLOCATABLE :: SEOF(:,:)
   REAL(DP),ALLOCATABLE :: SFSF(:,:)
   REAL(DP),ALLOCATABLE :: SFD(:)
   REAL(DP),ALLOCATABLE :: SFU(:,:)
   REAL(DP),ALLOCATABLE :: STVAR(:)
!   REAL(DP),ALLOCATABLE :: STSD(:)
   REAL(DP) :: SUM0   
   INTEGER RCODE
   INTEGER IDUMMY1
   REAL(DP),ALLOCATABLE :: RRKU2(:,:,:)
   REAL(DP),ALLOCATABLE :: RRKV2(:,:,:)
   REAL(DP),ALLOCATABLE :: RRKEL2(:,:)
   REAL(DP),ALLOCATABLE :: RRKT2(:,:,:)
   REAL(DP),ALLOCATABLE :: RRKS2(:,:,:)
!   REAL(DP),ALLOCATABLE :: STMEAN(:)

   TYPE(NCDIM), POINTER :: DIM
   TYPE(NCVAR), POINTER :: VAR

   INTEGER :: IRANGE(2),CNT, STATUS
   character(len=160) :: pathnfile 
   
   pathnfile = trim(INPUT_DIR)//trim(casename)//"_true_0001.nc"
   
   NCF => NEW_FILE()
   NCF%FNAME=trim(pathnfile)
   
   Call NC_OPEN(NCF)
   CALL NC_LOAD(NCF)

   ! GET THE TIME DIMENSION
   DIM => FIND_UNLIMITED(NCF,FOUND)
   IF(.not. FOUND) CALL FATAL_ERROR &
            & ("IN RRK_REF FILE",&
            & "FILE NAME: "//TRIM(NCF%FNAME),&
            &"COULD NOT FIND THE UNLIMITED DIMENSION")
   NTIMES = DIM%DIM
   

   ! TEST THE FILE START AND END
   TTEST = get_file_time(NCF,1)
   
   IF(REF_START_TIME < TTEST) CALL FATAL_ERROR &
        &("RRK_REF: REF_START_TIME IS BEFORE FIRST TIME IN FILE.",&
        & "FILE NAME:"//TRIM(NCF%FNAME))
   
   
   TTEST = get_file_time(NCF,NTIMES)
   IF(REF_END_TIME > TTEST) CALL FATAL_ERROR &
        &("RRK_REF: REF_START_TIME IS BEFORE FIRST TIME IN FILE.",&
        & "FILE NAME:"//TRIM(NCF%FNAME))


   ! GET THE START AND END INDEX

   IRANGE = 0

   CNT = 0
   DO WHILE (product(irange)==0)
      CNT = CNT +1

      IF (CNT > NTIMES) THEN
         CALL PRINT_TIME(REF_START_TIME,IPT,"RRK_REF START TIME")
         CALL PRINT_TIME(REF_END_TIME,IPT,"RRK_REF END TIME")
         

         CALL FATAL_ERROR&
           &("RRK_REF: CAN NOT FIND SPECIFIED START AND END TIME IN THE REFERENCE FILE!", &
           & "FILE NAME:"//TRIM(NCF%FNAME))
      END IF


      TTEST = get_file_time(NCF,CNT)
      IF (REF_START_TIME == TTEST) THEN
         IRANGE(1) = CNT
      END IF

      IF (REF_END_TIME == TTEST) THEN
         IRANGE(2) = CNT
      END IF

   END DO


   nrange = irange(2)-irange(1)+1

   IF(EL_ASSIM) THEN
      
      VAR => FIND_VAR(NCF,"zeta",FOUND)
      IF(.not. FOUND) CALL FATAL_ERROR &
           & ("IN RRK_REF FILE",&
            & "FILE NAME: "//TRIM(NCF%FNAME),&
            &"COULD NOT FIND THE VARIABLE 'zeta'")

      allocate(RRKEL_REF(MGL,nrange),stat=status)
      IF (STATUS /=0 ) CALL FATAL_ERROR("COULD NOT ALLOCATE MEMORY FOR RRKEL_REF")
      RRKEL_REF = 0.0_SP

      CALL nc_connect_avar(VAR,RRKEL_REF)

      CALL NC_READ_VAR(var,stkrng=irange,parallel=.false.)


   ENDIF
   
   IF(UV_ASSIM) THEN

      VAR => FIND_VAR(NCF,"u",FOUND)
      IF(.not. FOUND) CALL FATAL_ERROR &
           & ("IN RRK_REF FILE",&
            & "FILE NAME: "//TRIM(NCF%FNAME),&
            &"COULD NOT FIND THE VARIABLE 'u'")

      allocate(RRKU_REF(NGL,kbm1,nrange),stat=status)
      IF (STATUS /=0 ) CALL FATAL_ERROR("COULD NOT ALLOCATE MEMORY FOR RRKU_REF")
      RRKU_REF = 0.0_SP
      CALL nc_connect_avar(VAR,RRKU_REF)

      CALL NC_READ_VAR(var,stkrng=irange,parallel=.false.)

      VAR => FIND_VAR(NCF,"v",FOUND)
      IF(.not. FOUND) CALL FATAL_ERROR &
           & ("IN RRK_REF FILE",&
            & "FILE NAME: "//TRIM(NCF%FNAME),&
            &"COULD NOT FIND THE VARIABLE 'v'")

      allocate(RRKV_REF(NGL,kbm1,nrange),stat=status)
      IF (STATUS /=0 ) CALL FATAL_ERROR("COULD NOT ALLOCATE MEMORY FOR RRKV_REF")
      RRKV_REF = 0.0_SP
      CALL nc_connect_avar(VAR,RRKV_REF)

      CALL NC_READ_VAR(var,stkrng=irange,parallel=.false.)

      
   ENDIF
   
   IF(T_ASSIM) THEN
      VAR => FIND_VAR(NCF,"temp",FOUND)
      IF(.not. FOUND) CALL FATAL_ERROR &
           & ("IN RRK_REF FILE",&
            & "FILE NAME: "//TRIM(NCF%FNAME),&
            &"COULD NOT FIND THE VARIABLE 'temp'")

      allocate(RRKT_REF(MGL,kbm1,nrange),stat=status)
      IF (STATUS /=0 ) CALL FATAL_ERROR("COULD NOT ALLOCATE MEMORY FOR RRKT_REF")
      RRKT_REF = 0.0_SP
      CALL nc_connect_avar(VAR,RRKT_REF)

      CALL NC_READ_VAR(var,stkrng=irange,parallel=.false.)
      
   ENDIF
   
   IF(S_ASSIM) THEN
      VAR => FIND_VAR(NCF,"salinity",FOUND)
      IF(.not. FOUND) CALL FATAL_ERROR &
           & ("IN RRK_REF FILE",&
            & "FILE NAME: "//TRIM(NCF%FNAME),&
            &"COULD NOT FIND THE VARIABLE 'salinity'")

      allocate(RRKS_REF(MGL,kbm1,nrange),stat=status)
      IF (STATUS /=0 ) CALL FATAL_ERROR("COULD NOT ALLOCATE MEMORY FOR RRKS_REF")
      RRKS_REF = 0.0_SP
      CALL nc_connect_avar(VAR,RRKS_REF)

      CALL NC_READ_VAR(var,stkrng=irange,parallel=.false.)
      
   ENDIF
   
   CALL KILL_FILE(NCF)
   
   WRITE(IPT,*) 'Calculate the EOFs from the control run......'
   
   STDIM = 0
   IF(EL_ASSIM) STDIM = STDIM + MGL
   IF(UV_ASSIM) STDIM = STDIM + 2*NGL*KBM1
   IF(T_ASSIM)  STDIM = STDIM + MGL*KBM1
   IF(S_ASSIM)  STDIM = STDIM + MGL*KBM1
    
   SS_DIM = nrange
   LWORK4 = 5*SS_DIM
   LDVT   = 1   
   
   ALLOCATE(WORK4(LWORK4))             ;WORK4   = ZERO
   ALLOCATE(RKSF(STDIM,SS_DIM))        ;RKSF    = ZERO
   ALLOCATE(RKSF1(STDIM,SS_DIM))       ;RKSF1   = ZERO
   ALLOCATE(SEOF(STDIM,SS_DIM))        ;SEOF    = ZERO
   ALLOCATE(SFSF(SS_DIM,SS_DIM))       ;SFSF    = ZERO
   ALLOCATE(SFD(SS_DIM))               ;SFD     = ZERO
   ALLOCATE(SFU(SS_DIM,SS_DIM))        ;SFU     = ZERO
   ALLOCATE(STVAR(STDIM))              ;STVAR   = ZERO
   ALLOCATE(STSD(STDIM))               ;STSD    = ZERO
   ALLOCATE(RRKU(NGL,KBM1))            ;RRKU    = ZERO
   ALLOCATE(RRKV(NGL,KBM1))            ;RRKV    = ZERO
   ALLOCATE(RRKEL(MGL))                ;RRKEL   = ZERO
   ALLOCATE(RRKT(MGL,KBM1))            ;RRKT    = ZERO
   ALLOCATE(RRKS(MGL,KBM1))            ;RRKS    = ZERO
   ALLOCATE(RRKU2(NGL,KBM1,SS_DIM))    ;RRKU2   = ZERO
   ALLOCATE(RRKV2(NGL,KBM1,SS_DIM))    ;RRKV2   = ZERO
   ALLOCATE(RRKEL2(MGL,SS_DIM))        ;RRKEL2  = ZERO
   ALLOCATE(RRKT2(MGL,KBM1,SS_DIM))    ;RRKT2   = ZERO
   ALLOCATE(RRKS2(MGL,KBM1,SS_DIM))    ;RRKS2   = ZERO
   ALLOCATE(STMEAN(STDIM))             ;STMEAN  = ZERO  
   
   
   ! STORE IN RKSF

    
    DO I = 1, SS_DIM
  
      IDUMMY = 0      

      IF(EL_ASSIM) THEN      
        DO J = 1, MGL
          IDUMMY = IDUMMY + 1
	  RKSF(IDUMMY,I) =  RRKEL_REF(J,I)
        ENDDO 
      ENDIF

      IF(UV_ASSIM) THEN
        DO K = 1, KBM1
          DO J = 1, NGL
	     IDUMMY = IDUMMY + 1
	     RKSF(IDUMMY,I) = RRKU_REF(J,K,I)
	  ENDDO
        ENDDO	

        DO K = 1, KBM1
          DO J = 1, NGL
	     IDUMMY = IDUMMY + 1
	     RKSF(IDUMMY,I) = RRKV_REF(J,K,I)
	  ENDDO
        ENDDO     
      ENDIF

      IF(T_ASSIM) THEN
        DO K = 1, KBM1
          DO J = 1, MGL
	     IDUMMY = IDUMMY + 1
	     RKSF(IDUMMY,I) = RRKT_REF(J,K,I)
	  ENDDO
        ENDDO	
      ENDIF

      IF(S_ASSIM) THEN
	DO K = 1, KBM1
          DO J = 1, MGL
	     IDUMMY = IDUMMY + 1
	     RKSF(IDUMMY,I) = RRKS_REF(J,K,I)
	  ENDDO
        ENDDO	
      ENDIF
      
    ENDDO
    
      
! CALCULATE THE MEAN AND THE STANDARD DEVIATION

    DO I=1, STDIM
      SUM0=0.0_DP
      DO J=1, SS_DIM
        SUM0=SUM0+RKSF(I,J)
      ENDDO
     
      STMEAN(I) = SUM0/DBLE(SS_DIM)
      SUM0=0.0_DP
      DO J=1, SS_DIM
        RKSF1(I,J)=(RKSF(I,J)-STMEAN(I))
	SUM0=SUM0+RKSF1(I,J)**2.0_DP
      ENDDO
           
      STVAR(I)=SUM0
      STSD(I)=DSQRT(SUM0/DBLE(SS_DIM-1))
    ENDDO    
    

    IDUMMY  = 0
    IDUMMY1 = 0
    ELSD    = 0.0_DP
    USD     = 0.0_DP
    TSD     = 0.0_DP
    SSD     = 0.0_DP          

    IF(EL_ASSIM) THEN
      SUM0=0.0_DP
      DO I=1, MGL
        IDUMMY = IDUMMY + 1
        SUM0 = SUM0 + STVAR(IDUMMY)
      ENDDO
      
      ELSD=DSQRT(SUM0/DBLE(MGL)/DBLE(SS_DIM-1))

      DO I=1, MGL
        DO J=1, SS_DIM
          RKSF1(I,J)=RKSF1(I,J)/ELSD/DSQRT(DBLE(SS_DIM-1)) 
        ENDDO
      ENDDO
    
      IDUMMY1 = IDUMMY
    ENDIF

    IF(UV_ASSIM) THEN
      SUM0=0.0_DP
      DO K=1, 2*KBM1
        DO I=1, NGL
           IDUMMY = IDUMMY + 1
	   SUM0 = SUM0 +STVAR(IDUMMY)
        ENDDO
      ENDDO
      USD=DSQRT(SUM0/DBLE(NGL*KBM1)/DBLE(SS_DIM-1))
 
      DO I=IDUMMY1+1, IDUMMY1+2*NGL*KBM1 
        DO J=1, SS_DIM	    
	   RKSF1(I,J)=RKSF1(I,J)/USD/DSQRT(DBLE(KBM1))/DSQRT(DBLE(SS_DIM-1))         
        ENDDO 
      ENDDO
    
      IDUMMY1 = IDUMMY
    ENDIF  
      
    IF(T_ASSIM) THEN
      SUM0=0.0_DP
      DO K=1, KBM1
        DO I=1, MGL
           IDUMMY = IDUMMY + 1
	   SUM0 = SUM0 +STVAR(IDUMMY)
        ENDDO
      ENDDO
      TSD=DSQRT(SUM0/DBLE(MGL*KBM1)/DBLE(SS_DIM-1))
 
      DO I=IDUMMY1+1, IDUMMY1+MGL*KBM1 
        DO J=1, SS_DIM
	  RKSF1(I,J)=RKSF1(I,J)/TSD/DSQRT(DBLE(SS_DIM-1))
        ENDDO 
      ENDDO

      IDUMMY1 = IDUMMY
    ENDIF        
      
    IF(S_ASSIM) THEN
      SUM0=0.0_DP
      DO K=1, KBM1
        DO I=1, MGL
           IDUMMY = IDUMMY + 1
	   SUM0 = SUM0 +STVAR(IDUMMY)
        ENDDO
      ENDDO
      SSD=DSQRT(SUM0/DBLE(MGL*KBM1)/DBLE(SS_DIM-1))
 
      DO I=IDUMMY1+1, IDUMMY1+MGL*KBM1 
        DO J=1, SS_DIM
	   RKSF1(I,J)=RKSF1(I,J)/SSD/DSQRT(DBLE(SS_DIM-1))
        ENDDO 
      ENDDO

      IDUMMY1 = IDUMMY
    ENDIF       


! CALCULATE THE EOFs BY SVD OF RKSF1' RKSF1 = C LAMBDA C', RKSF1 RKSF1' = E LAMBDA E', E = RKSF1 C LAMBDA ^-1/2

    DO I=1, SS_DIM
      DO J=1, SS_DIM
        SUM0=0.0_DP
	DO K=1, STDIM
	  SUM0 = SUM0 + RKSF1(K,I)*RKSF1(K,J)
	ENDDO
	SFSF(I,J) =SUM0
      ENDDO
    ENDDO
 
    CALL DGESVD('A','N',SS_DIM,SS_DIM,SFSF,SS_DIM,SFD,SFU,SS_DIM,VT,LDVT,WORK4,LWORK4,RCODE)    
    
    OPEN(72,FILE=TRIM(OUTPUT_DIR)//'/rrktemp/'//'eigenvalue.dat')
    
    DO I=1, SS_DIM
      WRITE(72,'(I5,E15.7)') I, SFD(I)   
    ENDDO
    DO I=1, SS_DIM
      DO J=1, SS_DIM 
        SFU(I,J)=SFU(I,J)/DSQRT(SFD(J))
      ENDDO
    ENDDO

    DO I=1, STDIM
      DO J=1, SS_DIM
         SUM0=0.0_DP 
         DO K=1, SS_DIM
	    SUM0=SUM0+RKSF1(I,K)*SFU(K,J)
	 ENDDO
         SEOF(I,J)=SUM0 
      ENDDO
    ENDDO
    
    write(1000,*) (seof(i,1),i=1,stdim)
    
    DEALLOCATE(WORK4,RKSF,RKSF1,SFSF,SFD,SFU,STVAR)
    DEALLOCATE(RRKU2,RRKV2,RRKEL2,RRKT2,RRKS2)
    DEALLOCATE(RRKU,RRKV,RRKEL,RRKT,RRKS)

 END SUBROUTINE RRK_REF
 
!------------------------------------------------------------------------------|
!  Main program to calculate the reduced rank kalman gain matrix               |
!------------------------------------------------------------------------------|   
   
   SUBROUTINE RRK_RRK(CHOICE)

   USE LIMS
   USE CONTROL
   USE RRKVAL
   IMPLICIT NONE
    
    INTEGER SS_DIM
    INTEGER STDIM
    INTEGER I_EOF
    INTEGER I,J,K
    INTEGER CHOICE
!    INTEGER NLOC
    INTEGER RCODE
    INTEGER IT
    INTEGER ILAST
    INTEGER IDUMMY
    INTEGER SS_TOT
    INTEGER,ALLOCATABLE  :: STLOC(:)
    REAL(DP) ERR1
    REAL(DP),ALLOCATABLE :: KAL(:,:)
    REAL(DP),ALLOCATABLE :: HEOF(:,:)
    REAL(DP),ALLOCATABLE :: R(:,:)
    REAL(DP),ALLOCATABLE :: HTR(:,:)
    REAL(DP),ALLOCATABLE :: R2(:,:)
    REAL(DP),ALLOCATABLE :: PHT(:,:)
    REAL(DP),ALLOCATABLE :: RALPHA(:,:)
    REAL(DP),ALLOCATABLE :: BETA2(:,:)
    REAL(DP),ALLOCATABLE :: GAMMA(:,:)
    REAL(DP),ALLOCATABLE :: GAMMA2(:,:)
    REAL(DP),ALLOCATABLE :: TRANS2(:,:)
    
! WORK ARRAYS FOR LAPACK SUBROUTINES
    INTEGER LWORK, LWORK2, LDVT
    INTEGER,ALLOCATABLE :: IPIV(:)
    INTEGER,ALLOCATABLE :: IPIV2(:)
    REAL(DP),ALLOCATABLE :: WORK(:)
    REAL(DP),ALLOCATABLE :: WORK2(:)
    REAL(DP) :: VT
    

! CHARACTER STRINGS
    CHARACTER(LEN=80) KALNAME
    CHARACTER(LEN=80) FILENAME
    CHARACTER(LEN=4)  STRCYC
    CHARACTER(LEN=8)  CH8
    CHARACTER(LEN=80) INAME    
    CHARACTER(LEN=80) INAME2
    
!    REAL(DP),ALLOCATABLE :: SEOF(:,:)
!    REAL(DP),ALLOCATABLE :: STSD(:)
!    REAL(DP) :: ELSD, USD, TSD, SSD
    REAL(DP),ALLOCATABLE :: HU(:,:)
    REAL(DP),ALLOCATABLE :: HUL(:,:)
    REAL(DP),ALLOCATABLE :: EVAL(:)
    REAL(DP) :: RRKSUM
      
    STDIM = 0
    IF(EL_ASSIM) STDIM = STDIM + MGL
    IF(UV_ASSIM) STDIM = STDIM + 2*NGL*KBM1
    IF(T_ASSIM)  STDIM = STDIM + MGL*KBM1
    IF(S_ASSIM)  STDIM = STDIM + MGL*KBM1

    SS_DIM = RRK_NEOF
    LWORK  = 4*SS_DIM
    LWORK2 = 4*RRK_NOBSMAX
    SS_TOT = nrange
    LDVT   = 1 
    ILAST  = 10  ! THE NUMBER OF ITERATION IN DOUBLING ALGORITHM

! TEMPORARILY ALLOCATE ARRY TO ARRYS


   ALLOCATE(STLOC(RRK_NOBSMAX))            ;STLOC    = 0
   ALLOCATE(IPIV(SS_DIM))                  ;IPIV     = 0
   ALLOCATE(IPIV2(RRK_NOBSMAX))            ;IPIV2    = 0
   ALLOCATE(KAL(SS_DIM,RRK_NOBSMAX))       ;KAL      = ZERO
   ALLOCATE(HEOF(RRK_NOBSMAX,SS_DIM))      ;HEOF     = ZERO
   ALLOCATE(R(RRK_NOBSMAX,RRK_NOBSMAX))    ;R        = ZERO
   ALLOCATE(HTR(SS_DIM,RRK_NOBSMAX))       ;HTR      = ZERO
   ALLOCATE(R2(RRK_NOBSMAX,RRK_NOBSMAX))   ;R2       = ZERO
   ALLOCATE(PHT(SS_DIM,RRK_NOBSMAX))       ;PHT      = ZERO
   ALLOCATE(RALPHA(SS_DIM,SS_DIM))         ;RALPHA   = ZERO
   ALLOCATE(BETA2(SS_DIM,SS_DIM))          ;BETA2    = ZERO
   ALLOCATE(GAMMA(SS_DIM,SS_DIM))          ;GAMMA    = ZERO
   ALLOCATE(GAMMA2(SS_DIM,SS_DIM))         ;GAMMA2   = ZERO
   ALLOCATE(TRANS2(SS_DIM,SS_DIM))         ;TRANS2   = ZERO 
   ALLOCATE(WORK(LWORK))                   ;WORK     = ZERO
   ALLOCATE(WORK2(LWORK2))                 ;WORK2    = ZERO
!   ALLOCATE(SEOF(STDIM,SS_DIM))            ;SEOF     = ZERO
!   ALLOCATE(STSD(STDIM))                   ;STSD     = ZERO
   ALLOCATE(HU(RRK_NOBSMAX,SS_TOT))        ;HU       = ZERO      
   ALLOCATE(HUL(RRK_NOBSMAX,SS_TOT))       ;HUL      = ZERO
   ALLOCATE(EVAL(SS_TOT))                  ;EVAL     = ZERO


! END ALLOCATION

   IF(CHOICE == 1) THEN
     IF(MSR) WRITE(IPT,*) 'Starting perturb the base state in the eigenvector direction......'
     print *, 'call perturb' 
     CALL PERTURB

     DEALLOCATE(STLOC,IPIV,IPIV2,KAL,HEOF,R,HTR,R2,PHT,RALPHA,TRANS2)
     DEALLOCATE(BETA2,GAMMA,GAMMA2,WORK,WORK2,HU,HUL,EVAL) 

     RETURN
   ENDIF

   IF(CHOICE == 2) THEN
     IF(MSR) WRITE(IPT,*) 'Starting calculate the linearized model matrix in the reduced space......'  
! CALCULATE THE LINEARIZED MODEL MATRIX IN THE REDUCED SPACE
     CALL MREDUCED
     
     IF(MSR) THEN
! OUTPUT LINEAR TRANSITION MATRIX   
       FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Amatr.dat' 
       OPEN(INORRK,FILE=TRIM(FILENAME),FORM='UNFORMATTED')
       DO J=1, SS_DIM
         WRITE(INORRK) (TRANS(I,J), I=1, SS_DIM)    
       ENDDO 
       CLOSE(INORRK)

 ! ALSO IN ASCII
       FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Amatr1.txt'
       OPEN(INORRK, FILE=TRIM(FILENAME),FORM='FORMATTED')
       DO I=1, SS_DIM
         DO J=1, SS_DIM
           WRITE(INORRK,*) TRANS(J,I)
         ENDDO
       ENDDO
       CLOSE(INORRK)
     ENDIF  

     DEALLOCATE(STLOC,IPIV,IPIV2,KAL,HEOF,R,HTR,R2,PHT,RALPHA,TRANS2)
     DEALLOCATE(BETA2,GAMMA,GAMMA2,WORK,WORK2,HU,HUL,EVAL) 
!     DEALLOCATE(BETA2,GAMMA,GAMMA2,WORK,WORK2,SEOF,STSD,HU,HUL,EVAL)     
!     DEALLOCATE(STTEMP0,STTEMP1,STEOF,SDEOF,TRANS,RRKEL,RRKU,RRKV,RRKT,RRKS)
     RETURN
   ENDIF
       
   IF(CHOICE == 4) THEN
     CALL RRK_ALLOC_VAR
     WRITE(IPT,*) 'Calculate the Kalman gain matrix by doubling algorith......'
     KALNAME = TRIM(OUTPUT_DIR)//'/rrktemp/rrK.dat'

! READ THE EIGENVALUES AND SET THE INITIAL ERROR COVARIANCE   
     INAME = TRIM(OUTPUT_DIR)//'/rrktemp/eigenvalue.dat' 
     OPEN(INORRK,FILE=INAME,STATUS='OLD')
     
      WRITE(IPT,*) 'ss_tot:', ss_tot, ss_dim     
     
     DO I=1, SS_TOT
       READ(INORRK,*) J, EVAL(I)
     ENDDO
     
     DO I=1, SS_DIM
       DO J=1, SS_DIM
         GAMMA2(I,J) = 0._SP
       ENDDO
       GAMMA2(I,I) = EVAL(I)*RRK_PSCALE
     ENDDO
     CLOSE(INORRK)
   
! READ THE LINEAR TRANSITION MATRIX     
     FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Amatr.dat' 
     OPEN(INORRK,FILE=TRIM(FILENAME),FORM='UNFORMATTED')
     DO J=1, SS_DIM
       READ(INORRK) (TRANS(I,J), I=1, SS_DIM)    
     ENDDO 
     CLOSE(INORRK)
     
! READ THE OBSERVATION NUMBER AND LOCATION   
     CALL READOBS(STLOC,NLOC)
     
     print *, 'obs=', STLOC(1), NLOC
   
! SET RALPHA   
     DO I=1, SS_DIM
       DO J=1, SS_DIM
         TRANS2(I,J)=TRANS(J,I)
       ENDDO
!      WRITE(IPT,*) 'Amatr:', I, TRANS(I,I)     
     ENDDO

     print *, 'debug - 4'
       
! CALCULATE OBSERVATION MATRIX (EOF TO OBS TRANSFORMATION) H*D*EOFs
     DO J=1, SS_DIM
       DO I=1, STDIM  
          STTEMP1(I) = SEOF(I,J)
       ENDDO
       
       DO I=1, NLOC         
         HEOF(I,J)=STTEMP1(STLOC(I))*STSD(STLOC(I))
        WRITE(IPT,*) 'Heof:', I,J,HEOF(I,J)        
       ENDDO
     ENDDO
    

     FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/ObsOp.dat'
     OPEN(INORRK,FILE=FILENAME, FORM='UNFORMATTED')
     DO J=1, SS_DIM
       WRITE(INORRK) (HEOF(I,J), I=1, NLOC)    
     ENDDO 
     CLOSE(INORRK)     
     FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Heof.txt'
     OPEN(INORRK,FILE=FILENAME, FORM='FORMATTED')
     DO J=1, SS_DIM
       DO I=1, NLOC
         WRITE(INORRK,*) I,J,HEOF(I,J)
       ENDDO    
     ENDDO 
     CLOSE(INORRK)

! OUTPUT RALPHA MATRIX, WHICH IS TRANSPOSE OF TRANS
     FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Alpha.dat'
     OPEN(INORRK,FILE=FILENAME, FORM='UNFORMATTED')
     DO J=1, SS_DIM
       WRITE(INORRK) (TRANS2(I,J), I=1, SS_DIM)    
     ENDDO 
     CLOSE(INORRK)     

! OUTPUT GAMMA MATRIX
     FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Gamma.dat'
     OPEN(INORRK,FILE=FILENAME, FORM='UNFORMATTED')
     DO J=1, SS_DIM
       WRITE(INORRK) (GAMMA2(I,J), I=1, SS_DIM)    
     ENDDO 
     CLOSE(INORRK) 
     
! CALCULATE REPRESENTATIVE OBSERVATION ERROR DIRECTLY USING EOFs IN UNRESOLVED SUBSPACE
     DO I=SS_DIM+1, SS_TOT
!       CALL READEOF(I,2)
       DO J=1, STDIM  
          STTEMP1(J) = SEOF(J,I)
       ENDDO
       
       DO J=1, NLOC
          HU(J,I) = STTEMP1(STLOC(J))*STSD(STLOC(J))       
       ENDDO
     ENDDO
     DO I=1, NLOC
       RRKSUM=0.0_DP
       DO J=SS_DIM+1, SS_TOT
         HUL(I,J)=RRKSUM+HU(I,J)*EVAL(J)
       ENDDO
     ENDDO
     
     DO I=1, NLOC
       DO J=1, NLOC
          RRKSUM=0.0_DP
	  DO K=SS_DIM+1,SS_TOT
             RRKSUM = RRKSUM + HUL(I,K)*HU(J,K)
          ENDDO
          R(I,J) = RRKSUM
       ENDDO
       
! ASSUME THAT THE MEASUREMENT ERROR IS 1% OF THE ONE STANDARD DEVIATION OF THE CONTROL RUN
       IF(RRK_OPTION == 0) THEN
          R(I,I) = R(I,I) + (STSD(STLOC(I))*0.01_DP)**2.0_DP
       ELSE
          R(I,I) = R(I,I) + (STSD(STLOC(I))/DSQRT(DBLE(KBM1))*0.01_DP)**2.0_DP	  
       ENDIF
     ENDDO
   
     DO I=1, NLOC
       DO J=1, I-1
         R(I,J) = 0.5_DP*(R(I,J)+R(J,I))
         R(J,I) = R(I,J)
       ENDDO
     ENDDO
     FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Robs.dat'
     OPEN(INORRK,FILE=FILENAME, FORM='UNFORMATTED')
     DO J=1, NLOC
       WRITE(INORRK) (R(I,J), I=1, NLOC)    
     ENDDO 
     CLOSE(INORRK)  
     FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Robs.txt'
     OPEN(INORRK,FILE=FILENAME, FORM='FORMATTED')
     DO J=1, NLOC
       DO I=1, NLOC
         WRITE(INORRK,*) I,J, R(I,J)    
       ENDDO	 
     ENDDO 
     CLOSE(INORRK)   
     
! INVERT OBSERVATION ERROR COVARIANCE     
     CALL DGETRF(NLOC,NLOC,R,RRK_NOBSMAX,IPIV2,RCODE)
     IF(RCODE/=0) WRITE(IPT,*) 'error in computing LU factorization'     
     call DGETRI(NLOC,R,RRK_NOBSMAX,IPIV2,WORK2,LWORK2,RCODE)
     IF(RCODE/=0) WRITE(IPT,*) 'error in inverting the matrix'

! FORM BETA MATRIX = H_T*R**(-1)*H, WHERE H IS A NORMALIZED MATRIX H = H_ORIG*S
     CALL DGEMM('t','n',SS_DIM,NLOC,NLOC,1.0d0,HEOF,RRK_NOBSMAX,R,RRK_NOBSMAX,0.0d0,HTR,SS_DIM) 
     CALL DGEMM('n','n',SS_DIM,SS_DIM,NLOC,1.0d0,HTR,SS_DIM,HEOF,RRK_NOBSMAX,0.0d0,BETA2,SS_DIM)    
     FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Beta.dat'
     OPEN(INORRK,FILE=FILENAME, FORM='UNFORMATTED')
     DO J=1, SS_DIM
       WRITE(INORRK) (BETA2(I,J), I=1, SS_DIM)    
     ENDDO 
     CLOSE(INORRK)     
     
     WRITE(IPT,*) 'Running the doubling algorithm......'
   
     DO IT =1, ILAST 
       
! READ GAMMA FROM FILE
       FILENAME =  TRIM(OUTPUT_DIR)//'/rrktemp/Gamma.dat'    
       OPEN(INORRK,FILE=TRIM(FILENAME),FORM='UNFORMATTED')
       DO J=1, SS_DIM
         READ(INORRK) (GAMMA(I,J), I=1, SS_DIM)    
       ENDDO 
       CLOSE(INORRK)
                    
! READ BETA FROM FILE
       FILENAME =  TRIM(OUTPUT_DIR)//'/rrktemp/Beta.dat'    
       OPEN(INORRK,FILE=TRIM(FILENAME),FORM='UNFORMATTED')
       DO J=1, SS_DIM
         READ(INORRK) (BETA2(I,J), I=1, SS_DIM)    
       ENDDO 
       CLOSE(INORRK)
      
! COMPUTE EYE + BETA*GAMMA = STORE TEMPORARILY IN RALPHA (RALPHA = EYE + BETA*GAMMA)
         
       DO I=1, SS_DIM
         DO J=1, SS_DIM
           RALPHA(I,J)=0.0_DP
         ENDDO
         RALPHA(I,I)=1.0_DP
       ENDDO     
       CALL DGEMM('n','n',SS_DIM,SS_DIM,SS_DIM,1.0d0,BETA2,SS_DIM,GAMMA,SS_DIM,1.0d0,RALPHA,SS_DIM)
     
! COMPUTE INVERSE OF ABOBE RALPHA (=EYE+BETA*GAMMA) BY LAPACK ROUTINES AND STORE IT IN BETA (BETA=(EYE+BETA*GAMMA)**(-1))
       CALL DGETRF(SS_DIM,SS_DIM,RALPHA,SS_DIM,IPIV,RCODE)
  
       IF(RCODE/=0) WRITE(IPT,*) 'error in computing LU factorization'
       CALL DGETRI(SS_DIM,RALPHA,SS_DIM,IPIV,WORK,LWORK,RCODE)
       IF(RCODE/=0) WRITE(IPT,*) 'error in inverting the matrix'
       DO I=1, SS_DIM
         DO J=1, SS_DIM
           BETA2(I,J) = RALPHA(I,J)
	 ENDDO
       ENDDO
     
! READ RALPHA FROM FILE
       FILENAME =  TRIM(OUTPUT_DIR)//'/rrktemp/Alpha.dat'    
       OPEN(INORRK,FILE=TRIM(FILENAME),FORM='UNFORMATTED')
       DO J=1, SS_DIM
         READ(INORRK) (RALPHA(I,J), I=1, SS_DIM)    
       ENDDO 
       CLOSE(INORRK)

! COMPUTE PRODUCT (EYE+BETA*GAMMA)**(-1)*RALPHA (GAMMA = BETA'*RALPHA)
       CALL DGEMM('n','n',SS_DIM,SS_DIM,SS_DIM,1.0d0,BETA2,SS_DIM,RALPHA,SS_DIM,0.0d0,GAMMA,SS_DIM)  
          
! OUTPUT THIS PRODUCT TO TEMPORARY FILE
       FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/temp.02'
       OPEN(INORRK,FILE=TRIM(FILENAME),FORM='UNFORMATTED')
       DO J=1, SS_DIM
         WRITE(INORRK) (GAMMA(I,J), I=1, SS_DIM)    
       ENDDO 
       CLOSE(INORRK)     

! COMPUTE PRODUCT RALPHA*(EYE+BETA*GAMMA)**(-1) AND STORE IN GAMMA (GAMMA = RALPHA*BETA')
       CALL DGEMM('n','n',SS_DIM,SS_DIM,SS_DIM,1.0d0,RALPHA,SS_DIM,BETA2,SS_DIM,0.0d0,GAMMA,SS_DIM) 
   
! READ BACK OLD FILE
       FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Beta.dat'
       OPEN(INORRK,FILE=TRIM(FILENAME),FORM='UNFORMATTED')
       DO J=1, SS_DIM
         READ(INORRK) (Beta2(I,J), I=1, SS_DIM)    
       ENDDO 
       CLOSE(INORRK)  

! COMPUTE PRODUCT RALPHA*(EYE+BETA*GAMMA)**(-1)*BETA (RALPHA' = GAMMA'*BETA')
       CALL DGEMM('n','n',SS_DIM,SS_DIM,SS_DIM,1.0d0,GAMMA,SS_DIM,BETA2,SS_DIM,0.0d0,RALPHA,SS_DIM)

! READ BACK OLD RALPHA AND PUT IN GAMMA
       FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Alpha.dat'
       OPEN(INORRK,FILE=TRIM(FILENAME),FORM='UNFORMATTED')
       DO J=1, SS_DIM
         READ(INORRK) (GAMMA(I,J), I=1, SS_DIM)    
       ENDDO 
       CLOSE(INORRK)  

! COMPUTE BETA+RALPHA*(EYE+BETA*GAMMA)**(-1)*BETA*RALPHA_T (BETA = BETA+RALPHA'*GAMMA_T)
       CALL DGEMM('n','t',SS_DIM,SS_DIM,SS_DIM,1.0d0,RALPHA,SS_DIM,GAMMA,SS_DIM,1.0d0,BETA2,SS_DIM)

! MAKE SURE BETA IS SYMMETRIC
       DO I=1, SS_DIM
          DO J=1, I-1
             BETA2(I,J) = 0.5_DP*(BETA2(I,J)+BETA2(J,I))
             BETA2(J,I) = BETA2(I,J) 
	  ENDDO
       ENDDO

! SAVE THIS NEW BETA TO FILE
       FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Beta.dat'
       OPEN(INORRK,FILE=FILENAME, FORM='UNFORMATTED')
         DO J=1, SS_DIM
         WRITE(INORRK) (BETA2(I,J), I=1, SS_DIM)    
       ENDDO 
       CLOSE(INORRK)      

! READ BACK OLD GAMMA
       FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Gamma.dat'
       OPEN(INORRK,FILE=FILENAME, FORM='UNFORMATTED')
       DO J=1, SS_DIM
         READ(INORRK) (GAMMA(I,J), I=1, SS_DIM)    
       ENDDO 
       CLOSE(INORRK)     

! READ TEMPORARY FILE INTO RALPHA
       FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/temp.02'
       OPEN(INORRK,FILE=FILENAME, FORM='UNFORMATTED')
       DO J=1, SS_DIM
         READ(INORRK) (RALPHA(I,J), I=1, SS_DIM)    
       ENDDO 
       CLOSE(INORRK)       

! COMPUTE GAMMA*(EYE+BETA*GAMMA)**(-1)*RALPHA (BETA'=GAMMA*RALPHA')
       CALL DGEMM('n','n',SS_DIM,SS_DIM,SS_DIM,1.0d0,GAMMA,SS_DIM,RALPHA,SS_DIM,0.0d0,BETA2,SS_DIM)

! READ BACK OLD RALPHA
       FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Alpha.dat'
       OPEN(INORRK,FILE=FILENAME, FORM='UNFORMATTED')
       DO J=1, SS_DIM
         READ(INORRK) (RALPHA(I,J), I=1, SS_DIM)    
       ENDDO 
       CLOSE(INORRK)     

! COMPUTE GAMMA + RALPHA_T*GAMMA*(EYE+BETA*GAMMA)**(-1)*RALPHA (GAMMA = GAMMA+RALPHA_T*BETA')
       CALL DGEMM('t','n',SS_DIM,SS_DIM,SS_DIM,1.0d0,RALPHA,SS_DIM,BETA2,SS_DIM,1.0d0,GAMMA,SS_DIM)

! ENSURE GAMMA IS SYMMETRIC
       DO I=1, SS_DIM
          DO J=1, I-1
             GAMMA(I,J) = 0.5_DP*(GAMMA(I,J)+GAMMA(J,I))
             GAMMA(J,I) = GAMMA(I,J) 
	  ENDDO
       ENDDO

! SAVE THIS NEW GAMMA TO FILE
       FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Gamma.dat'
       OPEN(INORRK,FILE=FILENAME, FORM='UNFORMATTED')
       DO J=1, SS_DIM
         WRITE(INORRK) (GAMMA(I,J), I=1, SS_DIM)    
       ENDDO 
       CLOSE(INORRK)      
       FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Pfctr.txt'
       OPEN(INORRK,FILE=FILENAME, FORM='FORMATTED')
         DO J=1, SS_DIM
           DO I=1, SS_DIM
	      WRITE(INORRK,*) I,J,GAMMA(I,J)    
           ENDDO
	 ENDDO 
       CLOSE(INORRK)

! READ TEMPORARY FILE INTO GAMMA   
       FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/temp.02'
       OPEN(INORRK,FILE=FILENAME, FORM='UNFORMATTED')
       DO J=1, SS_DIM
         READ(INORRK) (GAMMA(I,J), I=1, SS_DIM)    
       ENDDO 
       CLOSE(INORRK)
   
! COMPUTE PRODUCT RALPHA*(EYE+BETA*GAMMA)**(-1)*RALPHA (BETA = RALPHA*GAMMA)  
       CALL DGEMM('n','n',SS_DIM,SS_DIM,SS_DIM,1.0d0,RALPHA,SS_DIM,GAMMA,SS_DIM,0.0d0,BETA2,SS_DIM)
   
! SAVE THIS NEW RALPHA TO FILE
       FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Alpha.dat'
       OPEN(INORRK,FILE=FILENAME, FORM='UNFORMATTED')
         DO J=1, SS_DIM
         WRITE(INORRK) (BETA2(I,J), I=1, SS_DIM)    
       ENDDO 
       CLOSE(INORRK)

!===================================================================================================
! END OF ITERATION FOR THE DOUBLING ALGORITHM          
     ENDDO   
!===================================================================================================

     WRITE(IPT,*) 'Setting up the Kalman gain matrix in the reduced space......'   

! READ IN PF FROM DOUBLING ALGORITHM 
     FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Gamma.dat'
     OPEN(INORRK,FILE=FILENAME, FORM='UNFORMATTED')
     DO J=1, SS_DIM
       READ(INORRK) (GAMMA(I,J), I=1, SS_DIM)    
     ENDDO 
     CLOSE(INORRK)
   
! READ IN OBSERVATION OPERATOR   
     FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/ObsOp.dat'
     OPEN(INORRK,FILE=FILENAME, FORM='UNFORMATTED')
     DO J=1, SS_DIM
       READ(INORRK) (HEOF(I,J), I=1, NLOC)    
     ENDDO 
     CLOSE(INORRK)
      
! READ OBSERVATION ERROR COVARIANCE      
     FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Robs.dat'
     OPEN(INORRK,FILE=FILENAME, FORM='UNFORMATTED')
     DO J=1, NLOC
       READ(INORRK) (R(I,J), I=1, NLOC)    
     ENDDO 
     CLOSE(INORRK)
   
! COMPUTE P*H^T (PHT = GAMMA * HEOF^T)   
     CALL DGEMM('n','t',SS_DIM,NLOC,SS_DIM,1.0d0,GAMMA,SS_DIM,HEOF,RRK_NOBSMAX,0.0d0,PHT,SS_DIM)

! COMPUTE H*P*H^T (R2 = H*P*HT)  
     CALL DGEMM('n','n',NLOC,NLOC,SS_DIM,1.0d0,HEOF,RRK_NOBSMAX,PHT,SS_DIM,0.0d0,R2,RRK_NOBSMAX)

! OUPUT BOTH FORCAST AND OBSERVATION ERROR STANDARD DEVIATION IN ASCII
     FILENAME =  TRIM(OUTPUT_DIR)//'/rrktemp/Fctobserr.txt'  
     OPEN(INORRK,FILE=FILENAME, FORM='FORMATTED')
     DO I=1, NLOC
       WRITE(INORRK,*) STLOC(I),DSQRT(R2(I,I)),DSQRT(R(I,I))
     ENDDO
     CLOSE(INORRK)
   
! ADD R + (H P H^t) ---> R   
     DO I = 1, NLOC
        DO J=1, NLOC
           R(J,I) = R(J,I) + R2(J,I)
        ENDDO
     ENDDO
   
! INVERT R (H*P*H^T+R)   
     CALL DGETRF(NLOC,NLOC,R,RRK_NOBSMAX,IPIV2,RCODE)
     IF(RCODE/=0) WRITE(IPT,*) 'error in computing LU factorization'
     CALL DGETRI(NLOC,R,RRK_NOBSMAX,IPIV2,WORK2,LWORK2,RCODE)
     IF(RCODE/=0) WRITE(IPT,*) 'error in inverting the matrix'
     
! COMPUTE KALMAN GAIN: K = P*HT*(H*P*H^T+R)^(-1)   
     CALL DGEMM('n','n',SS_DIM,NLOC,NLOC,1.0d0,PHT,SS_DIM,R,RRK_NOBSMAX,0.0d0,KAL,SS_DIM)
   
! OUTPUT RESULT   
     OPEN(INORRK,FILE=KALNAME, FORM='UNFORMATTED') 
     DO J=1, NLOC
       WRITE(INORRK) (KAL(I,J), I=1, SS_DIM)    
     ENDDO 
     CLOSE(INORRK)     
     FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/rrK.txt'
     OPEN(INORRK,FILE=FILENAME, FORM='FORMATTED') 
     DO J=1, NLOC
       DO I=1, SS_DIM
	 WRITE(INORRK,*) I,J,KAL(I,J)    
       ENDDO
     ENDDO 
     CLOSE(INORRK)
     
! OUTPUT KALMAN GAIN MATRIX IN THE FULL SPACE: D*ER*KAL
!     DO J=1, SS_DIM
!       CALL READEOF(J,2)
!       DO I=1, STDIM
!         SEOF(I,J)=STTEMP1(I)
!       ENDDO
!     ENDDO
     
     DO J=1, NLOC
       WRITE(STRCYC,'(I4.4)') J
       DO I=1, STDIM
          RRKSUM = 0.0_DP
          DO K=1, SS_DIM
             RRKSUM = RRKSUM + STSD(I)*SEOF(I,K)*KAL(K,J) 
          ENDDO 
	  STTEMP1(I) = RRKSUM
       ENDDO
       FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/'//'BigK'//STRCYC//'.txt'     
       OPEN(INORRK,FILE=FILENAME,FORM='FORMATTED')
       DO I=1, STDIM
          WRITE(INORRK,*) STTEMP1(I)
       ENDDO
       CLOSE(INORRK)
     ENDDO
     
! COMPUTE ANALYSIS ERROR COVARIANCE: P_AN1 = (I-K*H)*P (JUST AS DIAGNOSIS)
     CALL DGEMM('n','t',SS_DIM,SS_DIM,NLOC,-1.0d0,KAL,SS_DIM,PHT,SS_DIM,1.0d0,GAMMA,SS_DIM)
     FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Panlr.txt'
     OPEN(INORRK,FILE=FILENAME, FORM='FORMATTED') 
     DO J=1, SS_DIM
       DO I=1, SS_DIM
	 WRITE(INORRK,*) I,J,GAMMA(I,J)    
       ENDDO
     ENDDO 
     CLOSE(INORRK)
     
     CALL DGEMM('n','n',SS_DIM,SS_DIM,NLOC,1.0d0,KAL,SS_DIM,HEOF,RRK_NOBSMAX,0.0d0,RALPHA,SS_DIM)
     
! COMPUTE (I-KH)M AND ITS SINGULAR VALUE     
     DO I=1, SS_DIM
       DO J=1, SS_DIM
         IF(I/=J) THEN
            RALPHA(I,J) = (-1.0_DP)*RALPHA(I,J)
         ELSE
	    RALPHA(I,I) = 1.0_DP - RALPHA(I,I)
         ENDIF
       ENDDO
     ENDDO
     
     FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/I_KH.txt'
     OPEN(INORRK,FILE=FILENAME, FORM='FORMATTED') 
     DO J=1, SS_DIM
       DO I=1, SS_DIM
	 WRITE(INORRK,*) I,J,RALPHA(I,J)    
       ENDDO
     ENDDO 
     CLOSE(INORRK)     
     
     FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/Amatr.dat'
     OPEN(INORRK,FILE=FILENAME, FORM='UNFORMATTED')
     DO J=1, SS_DIM
       READ(INORRK) (TRANS(I,J), I=1, SS_DIM)    
     ENDDO 
     CLOSE(INORRK)     
     CALL DGEMM('n','n',SS_DIM,SS_DIM,SS_DIM,1.0d0,RALPHA,SS_DIM,TRANS,SS_DIM,0.0d0,RALPHA,SS_DIM)
    
     FILENAME = TRIM(OUTPUT_DIR)//'/rrktemp/M_KHM.txt'
     OPEN(INORRK,FILE=FILENAME, FORM='FORMATTED') 
     DO J=1, SS_DIM
       DO I=1, SS_DIM
	 WRITE(INORRK,*) I,J,RALPHA(I,J)    
       ENDDO
     ENDDO 
     CLOSE(INORRK) 
     
     DEALLOCATE(STLOC,IPIV,IPIV2,KAL,HEOF,R,HTR,R2,PHT,RALPHA,TRANS2)
!     DEALLOCATE(BETA2,GAMMA,GAMMA2,WORK,WORK2,SEOF,STSD,HU,HUL,EVAL)
     DEALLOCATE(BETA2,GAMMA,GAMMA2,WORK,WORK2,HU,HUL,EVAL)
!     DEALLOCATE(STTEMP0,STTEMP1,STEOF,SDEOF,TRANS,RRKU,RRKV,RRKEL,RRKT,RRKS)
   ENDIF
   
   RETURN
   END SUBROUTINE RRK_RRK


! UTILITIES PROGRAMS
 
  SUBROUTINE RRK_ALLOC_VAR

   USE CONTROL
   USE RRKVAL
   IMPLICIT NONE
   
   INTEGER STDIM
   INTEGER SS_DIM
   
   STDIM = 0
   IF(EL_ASSIM) STDIM = STDIM + MGL
   IF(UV_ASSIM) STDIM = STDIM + 2*NGL*KBM1
   IF(T_ASSIM)  STDIM = STDIM + MGL*KBM1
   IF(S_ASSIM)  STDIM = STDIM + MGL*KBM1
   
   SS_DIM = RRK_NEOF
   
! ALLOCATE ARRYS

   IF(.not.ALLOCATED(STTEMP0))  ALLOCATE(STTEMP0(STDIM))  ;STTEMP0   = ZERO
   IF(.not.ALLOCATED(STTEMP1))  ALLOCATE(STTEMP1(STDIM))  ;STTEMP1   = ZERO
   IF(.not.ALLOCATED(STEOF))    ALLOCATE(STEOF(STDIM))    ;STEOF     = ZERO
   IF(.not.ALLOCATED(SDEOF))    ALLOCATE(SDEOF(STDIM))    ;SDEOF     = ZERO
   IF(.not.ALLOCATED(TRANS))    ALLOCATE(TRANS(SS_DIM,SS_DIM)) ;TRANS     = ZERO
   IF(.not.ALLOCATED(RRKEL))    ALLOCATE(RRKEL(0:MGL))    ;RRKEL     = ZERO
   IF(.not.ALLOCATED(RRKU))     ALLOCATE(RRKU(0:NGL,KB))  ;RRKU      = ZERO
   IF(.not.ALLOCATED(RRKV))     ALLOCATE(RRKV(0:NGL,KB))  ;RRKV      = ZERO
   IF(.not.ALLOCATED(RRKT))     ALLOCATE(RRKT(0:MGL,KB))  ;RRKT      = ZERO
   IF(.not.ALLOCATED(RRKS))     ALLOCATE(RRKS(0:MGL,KB))  ;RRKS      = ZERO
   
! END ALLOCATION
   
   RETURN
  END SUBROUTINE RRK_ALLOC_VAR

!------------------------------------------------------------------------------|
!  PERTURB THE BASE STATE IN THE EIGENVECTOR DIRECTIONS                        |
!------------------------------------------------------------------------------|   
   
  SUBROUTINE PERTURB

#  if defined(WET_DRY)
     USE MOD_WD
#  endif
   USE LIMS
   USE ALL_VARS
   USE RRKVAL
   USE CONTROL
#  if defined (MULTIPROCESSOR)
     USE MOD_PAR
#  endif
   IMPLICIT NONE
    INTEGER IERR 
    INTEGER SS_DIM
    INTEGER STDIM
    INTEGER I_EOF
    INTEGER II,I,J
    INTEGER IDUMMY
    INTEGER I_START
    CHARACTER(LEN=80) IFILE
    CHARACTER(LEN=80) IFILE2
    CHARACTER(LEN=80) TEMP1FILE
    CHARACTER(LEN=80) TEMP2FILE
    CHARACTER(LEN=4)  FEOF
      
    STDIM = 0
    IF(EL_ASSIM) STDIM = STDIM + MGL
    IF(UV_ASSIM) STDIM = STDIM + 2*NGL*KBM1
    IF(T_ASSIM)  STDIM = STDIM + MGL*KBM1
    IF(S_ASSIM)  STDIM = STDIM + MGL*KBM1
    
    SS_DIM = RRK_NEOF

   ALLOCATE(STTEMP0(STDIM))                ;STTEMP0   = ZERO
   ALLOCATE(STTEMP1(STDIM))                ;STTEMP1   = ZERO
   ALLOCATE(STEOF(STDIM))                  ;STEOF     = ZERO
   ALLOCATE(SDEOF(STDIM))                  ;SDEOF     = ZERO
   ALLOCATE(TRANS(SS_DIM,SS_DIM))          ;TRANS     = ZERO        

    IDUMMY  = 0
    IF(EL_ASSIM) THEN
      DO I=1, MGL
        IDUMMY = IDUMMY + 1
        SDEOF(IDUMMY) = ELSD
      ENDDO
    ENDIF
    IF(UV_ASSIM) THEN
      DO I=1, 2*KBM1
        DO J=1, NGL 
          IDUMMY = IDUMMY + 1
          SDEOF(IDUMMY)=USD*DSQRT(DBLE(KBM1))     
        ENDDO
      ENDDO
    ENDIF
    IF(T_ASSIM) THEN
      DO I=1, KBM1
        DO J=1, MGL 
          IDUMMY = IDUMMY + 1
          SDEOF(IDUMMY)=TSD
        ENDDO
      ENDDO    
    ENDIF
    IF(S_ASSIM) THEN
      DO I=1, KBM1
        DO J=1, MGL 
          IDUMMY = IDUMMY + 1
          SDEOF(IDUMMY)=SSD
        ENDDO
      ENDDO    
    ENDIF
      print *, 'before rrk_rrk read restart, file name:',nc_start%fname ! should be marked
 
    print *, 'call readrestart' 
    CALL READRESTART(NC_START,REF_START_TIME)

!   IF(.not.ALLOCATED(RRKEL))    ALLOCATE(RRKEL3(0:MGL))    ;RRKEL3     = ZERO
!   IF(.not.ALLOCATED(RRKU))     ALLOCATE(RRKU3(0:NGL,KB))  ;RRKU3      = ZERO
!   IF(.not.ALLOCATED(RRKV3))     ALLOCATE(RRKV3(0:NGL,KB))  ;RRKV3      = ZERO
!   IF(.not.ALLOCATED(RRKT3))     ALLOCATE(RRKT3(0:MGL,KB))  ;RRKT3      = ZERO
!   IF(.not.ALLOCATED(RRKS3))     ALLOCATE(RRKS3(0:MGL,KB))  ;RRKS3      = ZERO
   ALLOCATE(RRKU(NGL,KB))            ;RRKU    = ZERO
   ALLOCATE(RRKV(NGL,KB))            ;RRKV    = ZERO
   ALLOCATE(RRKEL(MGL))                ;RRKEL   = ZERO
   ALLOCATE(RRKT(MGL,KB))            ;RRKT    = ZERO
   ALLOCATE(RRKS(MGL,KB))            ;RRKS    = ZERO

print *, allocated(rrku)

#  if defined (MULTIPROCESSOR)
if (par) then
CALL ACOLLECT(MYID,MSRID,NPROCS,NMAP,EL, RRKEL)
CALL ACOLLECT(MYID,MSRID,NPROCS,EMAP,U, RRKU)
CALL ACOLLECT(MYID,MSRID,NPROCS,EMAP,V, RRKV)
CALL ACOLLECT(MYID,MSRID,NPROCS,NMAP,T1, RRKT)
CALL ACOLLECT(MYID,MSRID,NPROCS,NMAP,S1, RRKS)
CALL MPI_BCAST(RRKEL,MGL,MPI_F,0,MPI_FVCOM_GROUP,IERR)
CALL MPI_BCAST(RRKU,NGL*KBM1,MPI_F,0,MPI_FVCOM_GROUP,IERR)
CALL MPI_BCAST(RRKV,NGL*KBM1,MPI_F,0,MPI_FVCOM_GROUP,IERR)
CALL MPI_BCAST(RRKT,MGL*KBM1,MPI_F,0,MPI_FVCOM_GROUP,IERR)
CALL MPI_BCAST(RRKS,MGL*KBM1,MPI_F,0,MPI_FVCOM_GROUP,IERR)

else
   RRKEL=EL
   RRKU=U
   RRKV=V
   RRKT=T
   RRKS=S
end if
#  elif
   RRKEL=EL
   RRKU=U
   RRKV=V
   RRKT=T
   RRKS=S
#  endif
    print *, 'call gr2st'
    CALL GR2ST(0)
    print *, 'finish gr2st'
    
! READ THE EOF AND PERTURB THE BASE STATE IN THIS DIRECTION
  if (msr) then
     print *, 'msr id' , myid
     open(343,file='seof',status='replace',form='unformatted')
     write(343) seof
     close(343)
  else
     print *, 'no msr id', myid
!    ALLOCATE(SEOF(STDIM,SS_DIM))   ! because the seof is calculated only the master cpu in rrk_ref
  end if
  
!write(myid+4000,*) myid 
# if defined (MULTIPROCESSOR)
IF(PAR)CALL MPI_BARRIER(MPI_FVCOM_GROUP,IERR)
     open(343,file='seof',status='old',form='unformatted',action='read')
     read(343) seof
     close(343)
IF(PAR)CALL MPI_BARRIER(MPI_FVCOM_GROUP,IERR)
! IF(PAR)CALL MPI_BCAST(SEOF,STDIM*SS_DIM,MPI_DOUBLE_precision,0,MPI_FVCOM_GROUP,IERR)
write(myid+3000,*) myid
# endif 
write(myid+2000,*) seof
   DO I_EOF=1, SS_DIM
!       CALL READEOF(I_EOF,1)
       DO I=1, STDIM  
          STEOF(I) = SEOF(I,I_EOF)
       ENDDO

! PERTURB THE BASE STATE IN THE DIRECTION OF THE I_EOF'TH EOF AND STORE IT IN THE MP1FILE (FOR RESTART)
       DO I=1, STDIM
          STTEMP1(I) = STTEMP0(I) + STEOF(I)*SDEOF(I)*DBLE(RRK_PSIZE)
           write(500+myid,*) STTEMP0(I), STEOF(I),SDEOF(I)
       ENDDO

     print *, 'call st2gr'
       CALL ST2GR
        do i=1,mgl
           write(400+myid,'(30f18.8)') (rrkt(i,j),j=1,kb)
        end do
      
   IF(SERIAL) THEN
     EL =  RRKEL 
     U  =  RRKU
     V  =  RRKV
     T1 =  RRKT
     S1 =  RRKS
   ELSE 
#  if defined (MULTIPROCESSOR)
     DO I=1, M
        EL(I)=RRKEL(NGID(I))
     ENDDO

     DO I=1, N
       DO J=1, KB
         U(I,J)=RRKU(EGID(I),J)
         V(I,J)=RRKV(EGID(I),J)
       ENDDO
     ENDDO

     DO I=1, M
       DO J=1, KB
         T1(I,J)=RRKT(NGID(I),J)
         S1(I,J)=RRKS(NGID(I),J)
       ENDDO
     ENDDO
#  endif
   ENDIF


! WRITE THE PERTURBED STATE IN TEMP1FILE
     CALL WRITERESTART(NCF_RRKRE,I_EOF)

   ENDDO


  RETURN
  END SUBROUTINE PERTURB

!------------------------------------------------------------------------------|
!  CALCULATE THE LINEARIZED MODEL IN THE REDUCED SPACE                         |
!------------------------------------------------------------------------------|   
   
  SUBROUTINE MREDUCED

   USE MOD_WD
   USE LIMS
   USE CONTROL
   USE RRKVAL
   IMPLICIT NONE
    
   INTEGER SS_DIM
   INTEGER STDIM
   INTEGER I_EOF
   INTEGER I,J
   INTEGER IDUMMY
   INTEGER I_START   
   CHARACTER(LEN=80) IFILE
   CHARACTER(LEN=80) IFILE2
   CHARACTER(LEN=80) TEMPFILE
   CHARACTER(LEN=4)  FEOF
   REAL(DP) ::  SUM0

   STDIM = 0
   IF(EL_ASSIM) STDIM = STDIM + MGL
   IF(UV_ASSIM) STDIM = STDIM + 2*NGL*KBM1
   IF(T_ASSIM)  STDIM = STDIM + MGL*KBM1
   IF(S_ASSIM)  STDIM = STDIM + MGL*KBM1

   SS_DIM = RRK_NEOF
   
   IDUMMY = 0
   IF(EL_ASSIM) THEN
     DO I=1, MGL
       IDUMMY = IDUMMY + 1
       SDEOF(IDUMMY) = ELSD
     ENDDO
   ENDIF
   IF(UV_ASSIM) THEN
     DO I =1, 2*KBM1
       DO J=1, NGL
         IDUMMY = IDUMMY + 1
         IF(RRK_OPTION == 0) THEN
            SDEOF(IDUMMY) = USD
         ELSE
            SDEOF(IDUMMY) = USD*DSQRT(DBLE(KBM1))
         ENDIF
       ENDDO
     ENDDO
   ENDIF
   IF(T_ASSIM) THEN
     DO I =1, KBM1
       DO J=1, MGL
         IDUMMY = IDUMMY + 1
         SDEOF(IDUMMY) = TSD
       ENDDO
     ENDDO
   ENDIF
   IF(S_ASSIM) THEN
     DO I =1, KBM1
       DO J=1, MGL
         IDUMMY = IDUMMY + 1
         SDEOF(IDUMMY) = SSD
       ENDDO
     ENDDO
   ENDIF   
   
   CALL READRESTART(NC_START,REF_START_TIME)
   CALL GR2ST(0)  
   
   
! READ THE FORECAST WHICH IS PROPAGATED FROM THE PERTURBED STATE  

   CALL NC_OPEN(NCF_RRKFCT)
 
   DO I_EOF =1, SS_DIM   
! READ THE EACH PERTURBED STATE AT THE END OF LINEARUZATION TIME STEP

   CALL READRESTART(NCF_RRKFCT,I_EOF) 

   CALL GR2ST(1)   
   
! PROJECT  EVOLVED PERTURBATION ONTO EOFs
     DO I=1, SS_DIM
       DO J=1, STDIM  
          STEOF(J) = SEOF(J,I)
       ENDDO
       SUM0 = 0.0_DP
       DO J=1, STDIM
         SUM0 = SUM0 + STEOF(J)/SDEOF(J)*(STTEMP1(J)-STTEMP0(J))
       ENDDO
       TRANS(I,I_EOF) = SUM0/DBLE(RRK_PSIZE)
           
     ENDDO
! EDN OF LOOP OVER EACH EOF   
   ENDDO
   
  END SUBROUTINE MREDUCED 


!=====================================================================================/
!  CONVERT THE STATE FROM THE GRID SPACE TO THE STATE VECTOR SPACE                    /
!=====================================================================================/
  SUBROUTINE GR2ST(OPT)
 
   USE LIMS
   USE ALL_VARS
   USE RRKVAL
   IMPLICIT NONE
   
    INTEGER IDUMMY, I, J
    INTEGER OPT

    IDUMMY=0
    
    IF(EL_ASSIM) THEN
      DO I=1, MGL
        IDUMMY = IDUMMY + 1
        IF (OPT == 0) THEN
          STTEMP0(IDUMMY) = RRKEL(I)
        ELSE
          STTEMP1(IDUMMY) = RRKEL(I)
        ENDIF	 	
      ENDDO
    ENDIF

    IF(UV_ASSIM) THEN
      DO I=1, KBM1
        DO J=1, NGL
          IDUMMY = IDUMMY + 1
          IF (OPT == 0) THEN
	    STTEMP0(IDUMMY) = RRKU(J,I)

	  ELSE
	    STTEMP1(IDUMMY) = RRKU(J,I)
	  ENDIF    
        ENDDO
      ENDDO
      DO I=1, KBM1
        DO J=1, NGL
          IDUMMY = IDUMMY + 1
          IF (OPT == 0) THEN
	    STTEMP0(IDUMMY) = RRKV(J,I)
	  ELSE
	    STTEMP1(IDUMMY) = RRKV(J,I)
	  ENDIF    
        ENDDO
      ENDDO		
    ENDIF  
  
    IF(T_ASSIM) THEN
      DO I=1, KBM1
        DO J=1, MGL
          IDUMMY = IDUMMY + 1
          IF (OPT == 0) THEN
	    STTEMP0(IDUMMY) = RRKT(J,I)
	  ELSE
	    STTEMP1(IDUMMY) = RRKT(J,I)
	  ENDIF    
        ENDDO
      ENDDO
    ENDIF
  
    IF(S_ASSIM) THEN
      DO I=1, KBM1
        DO J=1, MGL
          IDUMMY = IDUMMY + 1
          IF (OPT == 0) THEN
	    STTEMP0(IDUMMY) = RRKS(J,I)
	  ELSE
	    STTEMP1(IDUMMY) = RRKS(J,I)
	  ENDIF    
        ENDDO
      ENDDO
    ENDIF

  RETURN
  END SUBROUTINE GR2ST

!=====================================================================================/
!  CONVERT THE STATE FROM THE STATE VECTOR SPACE TO THE GRID SPACE                    /
!=====================================================================================/
  SUBROUTINE ST2GR
 
   USE LIMS
   USE ALL_VARS
   USE RRKVAL
   IMPLICIT NONE
    
    INTEGER IDUMMY, I, J

    IDUMMY=0
    IF(EL_ASSIM) THEN
      DO I=1, MGL
        IDUMMY = IDUMMY + 1
        RRKEL(I) = STTEMP1(IDUMMY) 
      ENDDO
    ENDIF
    IF(UV_ASSIM) THEN
      DO I=1, KBM1
        DO J=1, NGL
          IDUMMY = IDUMMY + 1
          RRKU(J,I) = STTEMP1(IDUMMY) 
        ENDDO
      ENDDO
      DO I=1, KBM1
        DO J=1, NGL
          IDUMMY = IDUMMY + 1
          RRKV(J,I) = STTEMP1(IDUMMY)
        ENDDO
      ENDDO
    ENDIF
    IF(T_ASSIM) THEN
      DO I=1, KBM1
        DO J=1, MGL
          IDUMMY = IDUMMY + 1
          RRKT(J,I) = STTEMP1(IDUMMY) 
        ENDDO
      ENDDO
    ENDIF
    IF(S_ASSIM) THEN
      DO I=1, KBM1
        DO J=1, MGL
          IDUMMY = IDUMMY + 1
          RRKS(J,I) = STTEMP1(IDUMMY) 
        ENDDO
      ENDDO
    ENDIF
    
  RETURN
  END SUBROUTINE ST2GR


!=====================================================================================/
!  DETERMINE IF NODES/ELEMENTS ARE WET OR DRY                                         /
!=====================================================================================/
  SUBROUTINE WET_JUDGE_EL

   USE MOD_PREC
   USE ALL_VARS
#  if defined (MULTIPROCESSOR)
   USE MOD_PAR
#  endif
#  if defined(WET_DRY)
   USE MOD_WD
#  endif
   IMPLICIT NONE
   REAL(SP) :: DTMP
   INTEGER  :: ITA_TEMP
   INTEGER  :: I,IL,IA,IB,K1,K2,K3,K4,K5,K6

#  if defined(WET_DRY)

!
!--Determine If Node Points Are Wet/Dry Based on Depth Threshold---------------!
!
   ISWETN = 1
   DO I = 1, M
     DTMP = H(I) + EL(I)
     IF((DTMP - MIN_DEPTH) < 1.0E-5_SP) ISWETN(I) = 0
   END DO

!
!--Determine if Cells are Wet/Dry Based on Depth Threshold---------------------!
!
   ISWETC = 1
   DO I = 1, N
     DTMP =  MAX(EL(NV(I,1)),EL(NV(I,2)),EL(NV(I,3)))  + &
             MIN(  H(NV(I,1)),  H(NV(I,2)),  H(NV(I,3)))
     IF((DTMP - MIN_DEPTH) < 1.0E-5_SP) ISWETC(I) = 0
   END DO

!
!--A Secondary Condition for Nodal Dryness-(All Elements Around Node Are Dry)--!
!
   DO I = 1, M
     IF(SUM(ISWETC(NBVE(I,1:NTVE(I)))) == 0)  ISWETN(I) = 0
   END DO

!
!--Adjust Water Surface So It Does Not Go Below Minimum Depth------------------!
!
   EL = MAX(EL,-H + MIN_DEPTH)

!
!--Recompute Element Based Depths----------------------------------------------!
!
   DO I = 1, N
     EL1(I) = ONE_THIRD*(EL(NV(I,1))+EL(NV(I,2))+EL(NV(I,3)))
   END DO

!
!--Extend Element/Node Based Wet/Dry Flags to Domain Halo----------------------!
!
#  if defined (MULTIPROCESSOR)
   IF(PAR)THEN

     CALL AEXCHANGE(EC,MYID,NPROCS,ISWETC)
     CALL AEXCHANGE(NC,MYID,NPROCS,ISWETN)

   END IF

#  endif
#  endif

   RETURN

  END SUBROUTINE WET_JUDGE_EL

!!$!==============================================================================|
!!$!   DUMP WET/DRY FLAG DATA FOR RESTART                                         |
!!$!==============================================================================|
!!$
!!$   SUBROUTINE WD_DUMP_EL(INF,I_START)
!!$
!!$!------------------------------------------------------------------------------|
!!$
!!$   USE ALL_VARS
!!$#  if defined (MULTIPROCESSOR)
!!$   USE MOD_PAR
!!$#  endif
!!$#  if defined(WET_DRY)
!!$   USE MOD_WD
!!$#  endif
!!$
!!$   IMPLICIT NONE
!!$   INTEGER, ALLOCATABLE,DIMENSION(:) :: NTEMP1,NTEMP2
!!$   INTEGER I, INF
!!$   INTEGER I_START
!!$!==============================================================================|
!!$
!!$#  if defined(WET_DRY)
!!$
!!$   IF(MSR)THEN
!!$     REWIND(INF)
!!$     WRITE(INF,*) I_START
!!$     WRITE(INF,*) NGL,MGL
!!$   END IF
!!$
!!$   IF(SERIAL)THEN
!!$     WRITE(INF,*) (ISWETC(I), I=1,N)
!!$     WRITE(INF,*) (ISWETN(I), I=1,M)
!!$   ELSE
!!$   ALLOCATE(NTEMP1(NGL),NTEMP2(MGL))
!!$#  if defined (MULTIPROCESSOR)
!!$   CALL IGATHER(LBOUND(ISWETC,1),UBOUND(ISWETC,1),N,NGL,1,MYID,NPROCS,EMAP,ISWETC,NTEMP1)
!!$   CALL IGATHER(LBOUND(ISWETN,1),UBOUND(ISWETN,1),M,MGL,1,MYID,NPROCS,NMAP,ISWETN,NTEMP2)
!!$   IF(MSR)THEN
!!$     WRITE(INF,*) (NTEMP1(I), I=1,NGL)
!!$     WRITE(INF,*) (NTEMP2(I), I=1,MGL)
!!$   END IF
!!$   DEALLOCATE(NTEMP1,NTEMP2)
!!$#  endif
!!$   END IF
!!$
!!$   CLOSE(INF)
!!$
!!$#  endif
!!$
!!$   RETURN
!!$   END SUBROUTINE WD_DUMP_EL

  SUBROUTINE READOBS(STLOC,NLOC)
    
   USE LIMS
   USE CONTROL
   IMPLICIT NONE
   
   INTEGER ::  NUM  = 0 
   INTEGER ::  NLOC
   INTEGER ::  SWITCH = 0
   INTEGER ::  J,K
   INTEGER ::  IDUMMY
   INTEGER ::  TMP
   INTEGER ::  ioerr
   CHARACTER(LEN=80) FILENAME
   CHARACTER(LEN=24) HEADINFO
   INTEGER STLOC(RRK_NOBSMAX)
   INTEGER LAY(RRK_NOBSMAX)

   FILENAME = TRIM(INPUT_DIR)//"/"//trim(casename)//"_assim_rrkf.nml"
    
   OPEN(73,FILE=TRIM(FILENAME),FORM='FORMATTED')
   REWIND(73)

   NLOC = 0
   IDUMMY = 0
   STLOC = ZERO

   DO WHILE(.TRUE.)
     READ(73,'(A24)',IOSTAT=ioerr) HEADINFO
     IF(ioerr .ne. 0) EXIT
     IF(SWITCH/=1) THEN
       IF(HEADINFO=='!=== READ IN OBSERVATION') SWITCH = 1
       CYCLE
     ENDIF
   
     IF(TRIM(HEADINFO)=='!EL') THEN
       IF(EL_OBS) THEN
         READ(73,*) NUM
         NLOC = NLOC + NUM
         IF(NLOC>RRK_NOBSMAX) THEN
           WRITE(IPT,*) 'not enough storage for observations:', 'Nloc=', Nloc, 'Nobsmax=', RRK_NOBSMAX
           CALL PSTOP
         ENDIF
         READ(73,*)  (STLOC(K), K=1,NLOC)	
       ENDIF  
     
       IF(EL_ASSIM) THEN
         IDUMMY = IDUMMY + MGL
       ENDIF
     ENDIF
     
     IF(TRIM(HEADINFO)=='!UV') THEN
       IF(UV_OBS) THEN       
         READ(73,*) NUM
         NLOC = NLOC + NUM
         IF(NLOC+NUM>RRK_NOBSMAX) THEN
           WRITE(IPT,*) 'not enough storage for observations:', 'Nloc=', Nloc+num, 'Nobsmax=', RRK_NOBSMAX
           CALL PSTOP
         ENDIF
         READ(73,*)  (STLOC(K), K=NLOC-NUM+1,NLOC)
         READ(73,*)  (LAY(K),   K=NLOC-NUM+1,NLOC)
         DO K=NLOC-NUM+1, NLOC
	         STLOC(K)=STLOC(K)+IDUMMY+NGL*(LAY(K)-1)
         ENDDO   
         IDUMMY = IDUMMY + NGL*KBM1
	 
         NLOC = NLOC + NUM
         DO K=NLOC-NUM+1, NLOC
           STLOC(K)=STLOC(K-NUM)+NGL*KBM1+NGL*(LAY(K-NUM)-1)	   
         ENDDO
       ENDIF
       
       IF(UV_ASSIM) THEN  
         IDUMMY = IDUMMY + NGL*KBM1
       ENDIF
     ENDIF
     
     IF(TRIM(HEADINFO)=='!T') THEN
       IF(T_OBS) THEN
         READ(73,*) NUM
         NLOC = NLOC + NUM
         IF(NLOC>RRK_NOBSMAX) THEN
           WRITE(IPT,*) 'not enough storage for observations:', 'Nloc=', Nloc, 'Nobsmax=', RRK_NOBSMAX
           CALL PSTOP
         ENDIF
         READ(73,*)  (STLOC(K), K=NLOC-NUM+1,NLOC)       
         READ(73,*)  (LAY(K),   K=NLOC-NUM+1,NLOC)      
         DO K=NLOC-NUM+1, NLOC
           STLOC(K)=STLOC(K)+IDUMMY+MGL*(LAY(K)-1)
         ENDDO     
       ENDIF   
       
       IF(T_ASSIM) THEN
         IDUMMY = IDUMMY + MGL*KBM1
       ENDIF
     ENDIF
         
     IF(TRIM(HEADINFO)=='!S') THEN
       IF(S_OBS) THEN
         READ(73,*) NUM
         NLOC = NLOC + NUM
         IF(NLOC>RRK_NOBSMAX) THEN
           WRITE(IPT,*) 'not enough storage for observations:', 'Nloc=', Nloc, 'Nobsmax=', RRK_NOBSMAX
           CALL PSTOP
         ENDIF
         READ(73,*)  (STLOC(K),K=NLOC-NUM+1,NLOC)        
         READ(73,*)  (LAY(K),  K=NLOC-NUM+1,NLOC)    
         DO K=NLOC-NUM+1, NLOC
           STLOC(K)=STLOC(K)+IDUMMY+MGL*(LAY(K)-1)
         ENDDO   
       ENDIF
       
       IF(S_ASSIM) THEN
         IDUMMY = IDUMMY + MGL*KBM1 
       ENDIF 
     ENDIF
   ENDDO     
   CLOSE(73)
    
   RETURN 
  END SUBROUTINE READOBS  

  SUBROUTINE READRESTART_STATE(NCF,STATE)
    USE MOD_INPUT
    USE MOD_STARTUP
    IMPLICIT NONE
    TYPE(NCFILE), POINTER :: NCF, NCF_tmP
    INTEGER :: STATE
    NCF_tmP=>NC_START
    IF(.not. associated(NCF)) CALL FATAL_ERROR&
         & ("MOD_RRK: READRESTART: THE FILE POINTER PASSED IS NOT ASSOCIATED?")  
    NCF%FTIME%PREV_STKCNT = STATE
    NC_START=>NCF
    
    

    CALL READ_SSH

    IF(WETTING_DRYING_ON) CALL READ_WETDRY

    CALL READ_UV

    CALL READ_TURB

    CALL READ_TS

    NC_START=>NCF_tmP

  END SUBROUTINE READRESTART_STATE

  SUBROUTINE READRESTART_TIME(NCF,NOW)
    USE MOD_INPUT
    USE MOD_STARTUP
    IMPLICIT NONE
    TYPE(NCFILE), POINTER :: NCF
    TYPE(TIME) :: NOW,TTEST
    INTEGER :: STATE,NSTATE
    TYPE(NCDIM), POINTER :: DIM
    LOGICAL FOUND
    
    IF(.not. associated(NCF)) CALL FATAL_ERROR&
         & ("MOD_RRK: READRESTART: THE FILE POINTER PASSED IS NOT ASSOCIATED?")

    ! GET THE TIME DIMENSION
   DIM => FIND_UNLIMITED(NCF,FOUND)
   IF(.not. FOUND) CALL FATAL_ERROR &
            & ("IN RESTART FILE",&
            & "FILE NAME: "//TRIM(NCF%FNAME),&
            &"COULD NOT FIND THE UNLIMITED DIMENSION")
   NSTATE = DIM%DIM
   FOUND = .false.
   DO STATE=1,NSTATE
      
      TTEST = get_file_time(NCF,STATE)
      IF(TTEST == NOW) THEN
         FOUND = .true.
         exit
      END IF
      
   END DO
   
   IF(FOUND)THEN
      CALL READRESTART_STATE(NCF,STATE)
   ELSE
      CALL PRINT_TIME(NOW,IPT,"TIME REQUSTED")
      CALL FATAL_ERROR&
           &("COULD NOT FIND CORRECT TIME IN RESTART FILE",&
           &"FILE NAME: "//TRIM(NCF%FNAME))
 
  END IF

  END SUBROUTINE READRESTART_TIME

  SUBROUTINE WRITERESTART(NCF,STATE)
    IMPLICIT NONE
    TYPE(NCFILE), POINTER :: NCF
    INTEGER :: STATE

    IF(.not. associated(NCF)) CALL FATAL_ERROR&
         & ("MOD_RRK: WRITERESTART: THE FILE POINTER PASSED IS NOT ASSOCIATED?")

    NCF%FTIME%NEXT_STKCNT = state
    CALL NC_WRITE_FILE(NCF,LOCAL_DISK)
    
  END SUBROUTINE WRITERESTART


!======================================================|
   FUNCTION SETUP_RESTART(FNAME) RESULT(NCF)
    USE MOD_NCDIO
    IMPLICIT NONE
    CHARACTER(LEN=*) :: FNAME    
    TYPE(GRID), SAVE :: MYGRID
    TYPE(NCFILE), POINTER :: NCF
    TYPE(NCFILE), POINTER :: NCF2
    IF(DBG_SET(DBG_SBR)) WRITE(IPT,*) "START SETUP_RSTFILE"

    IF(DBG_SET(DBG_LOG)) THEN
       
       write(IPT,*)"!--------------------------------------------------"
       write(IPT,*)"! SETTING UP RESTART FILE OUTPUT..."
   END IF
    NCF => NEW_FILE()
    NCF%FNAME=trim(fname)
    print *, 'running setup_restart:',NCF%FNAME ! should be mark
    NCF%FTIME=>NEW_FTIME()

    CALL SET_FVCOM_GRID(MYGRID)

    CALL DEFINE_DIMENSIONS(MYGRID)

    NCF2 => GRID_FILE_OBJECT(MYGRID)
    NCF => ADD(NCF,NCF2)
!!$    NCF => ADD(NCF,GRID_FILE_OBJECT(MYGRID) )

    NCF2 => TIME_FILE_OBJECT()
    NCF => ADD(NCF,NCF2)
!!$    NCF => ADD(NCF,TIME_FILE_OBJECT() )

    NCF2 => ZETA_FILE_OBJECT()
    NCF => ADD(NCF,NCF2)
!!$    NCF => ADD(NCF,ZETA_FILE_OBJECT() )

    NCF2 => FILE_DATE_OBJECT()
    NCF => ADD(NCF,NCF2)
!!$    NCF => ADD(NCF,FILE_DATE_OBJECT() )

    NCF2 => VELOCITY_FILE_OBJECT()
    NCF => ADD(NCF,NCF2)
!!$    NCF => ADD(NCF,VELOCITY_FILE_OBJECT() )

    NCF2 => AVERAGE_VEL_FILE_OBJECT()
    NCF => ADD(NCF,NCF2)
!!$    NCF => ADD(NCF,AVERAGE_VEL_FILE_OBJECT() )

    NCF2 => VERTICAL_VEL_FILE_OBJECT()
    NCF => ADD(NCF,NCF2)
!!$    NCF => ADD(NCF,VERTICAL_VEL_FILE_OBJECT() )

    NCF2 => TURBULENCE_FILE_OBJECT()
    NCF => ADD(NCF,NCF2)
!!$    NCF => ADD(NCF,TURBULENCE_FILE_OBJECT() )
    
    NCF2 => SALT_TEMP_FILE_OBJECT()
    NCF => ADD(NCF,NCF2)
!!$    NCF => ADD(NCF,SALT_TEMP_FILE_OBJECT() )

    NCF2 => DENSITY_FILE_OBJECT()
    NCF => ADD(NCF,NCF2)
!!$    NCF => ADD(NCF,DENSITY_FILE_OBJECT() )

    NCF2 => RESTART_EXTRAS_FILE_OBJECT()
    NCF => ADD(NCF,NCF2)
!!$    NCF => ADD(NCF,RESTART_EXTRAS_FILE_OBJECT() )
    
    IF(WETTING_DRYING_ON) THEN
       NCF2 => WET_DRY_FILE_OBJECT()
       NCF => ADD(NCF, NCF2)
!!$       NCF => ADD(NCF, WET_DRY_FILE_OBJECT() )
    END IF

#   if defined (BioGen)
    NCF2 => BIO_FILE_OBJECT()
    NCF => ADD(NCF,NCF2)
!!$    NCF => ADD(NCF,BIO_FILE_OBJECT() )
#   endif      

#   if defined (WATER_QUALITY)
    NCF2 => WQM_FILE_OBJECT()
    NCF => ADD(NCF,NCF2)
!!$    NCF => ADD(NCF,WQM_FILE_OBJECT() )
#   endif      

#   if defined (SEDIMENT)
    IF(SEDIMENT_MODEL)THEN
      NCF2 => SED_FILE_OBJECT()
      NCF => ADD(NCF,NCF2)
!!$      NCF => ADD(NCF,SED_FILE_OBJECT() )
    ENDIF
#   endif      

#   if defined (WAVE_CURRENT_INTERACTION)      
    NCF2 => WAVE_PARA_FILE_OBJECT()
    NCF => ADD(NCF,NCF2) 
!!$    NCF => ADD(NCF,WAVE_PARA_FILE_OBJECT() ) 

    NCF2 => WAVE_STRESS_FILE_OBJECT()
    NCF => ADD(NCF,NCF2) 
!!$    NCF => ADD(NCF,WAVE_STRESS_FILE_OBJECT() ) 
#   endif  


#   if defined (NH)
    NCF2 => NH_RST_FILE_OBJECT()
    NCF => ADD(NCF,NCF2)
!!$    NCF => ADD(NCF,NH_RST_FILE_OBJECT() )
#   endif

# if defined (ICE)
      !-----------------------------------------------------------------
      ! state variables
    NCF2 => ICE_RESTART_STATE_FILE_OBJECT()
    NCF => ADD(NCF,NCF2)
!!$    NCF => ADD(NCF,ICE_RESTART_STATE_FILE_OBJECT() )
      !-----------------------------------------------------------------
      ! velocity
      !-----------------------------------------------------------------
    NCF2 => ICE_VEL_FILE_OBJECT()
    NCF => ADD(NCF,NCF2)
!!$    NCF => ADD(NCF,ICE_VEL_FILE_OBJECT() )
      !-----------------------------------------------------------------
      ! fresh water, salt, and heat flux
      !-----------------------------------------------------------------
    NCF2 => ICE_EXTRA_FILE_OBJECT()
    NCF => ADD(NCF,NCF2)
!!$    NCF => ADD(NCF,ICE_EXTRA_FILE_OBJECT() )
# endif

    IF (STARTUP_TYPE /= "crashrestart") THEN
       CALL NC_WRITE_FILE(NCF)
       NCF%FTIME%NEXT_STKCNT = 1
    ELSE
       NCF%CONNECTED = .TRUE.
    END IF

    CALL KILL_DIMENSIONS

    IF(DBG_SET(DBG_SBR)) WRITE(IPT,*) "END SETUP_RSTFILE"
  END FUNCTION SETUP_RESTART




# endif
END MODULE MOD_RRK

